Conditional Spatial Plotting with GRIB2 data

[edrewitz]

Hi all!

I am trying to plot on a cartopy map in one shaded color where the following is true: Temperature >=75F, Relative Humidity <=25% and Wind Speed >=15mph. Does anyone have any tips on how I can compare all those 3 parameters and plot them as one shaded area where those 3 conditions hold true? To be a little more specific, I am plotting the current RTMA data and I used MetPy to calculate RH from temperature and dewpoint. Thanks!

Here is the section code I have so far but it doesn't work at the end where I am trying to use the where() method in xarray:

P.S. I set the temperature to 15F for now even though that isn't RFW criteria since I am testing this out for Alaska where there won't be any RFW criteria for several months :)

        def get_RTMA_Data_red_flag_warning(current_time):

            r'''

            THIS FUNCTION RETRIEVES THE RTMA DATA TO DETERMINE CURRENT AREAS MEETING RED FLAG WARNING CRITERIA

            (C) METEOROLOGIST ERIC J. DREWITZ 2024

            '''
            
            times = []
            for i in range(0, 5):
                time = pd.to_datetime(current_time - timedelta(hours=i))
                times.append(time)

            url_0 = 'http://nomads.ncep.noaa.gov:80/dods/akrtma/akrtma'+times[0].strftime('%Y%m%d')+'/akrtma_anl_'+times[0].strftime('%H')+'z'
            url_1 = 'http://nomads.ncep.noaa.gov:80/dods/akrtma/akrtma'+times[1].strftime('%Y%m%d')+'/akrtma_anl_'+times[1].strftime('%H')+'z'
            url_2 = 'http://nomads.ncep.noaa.gov:80/dods/akrtma/akrtma'+times[2].strftime('%Y%m%d')+'/akrtma_anl_'+times[2].strftime('%H')+'z'
            url_3 = 'http://nomads.ncep.noaa.gov:80/dods/akrtma/akrtma'+times[3].strftime('%Y%m%d')+'/akrtma_anl_'+times[3].strftime('%H')+'z'
            url_4 = 'http://nomads.ncep.noaa.gov:80/dods/akrtma/akrtma'+times[4].strftime('%Y%m%d')+'/akrtma_anl_'+times[4].strftime('%H')+'z'

            try:
                ds = xr.open_dataset(url_0, engine='netcdf4')
                print("Data was successfully retrieved for " + times[0].strftime('%m/%d/%Y %HZ'))
                strtime = times[0]
                
            except Exception as a:
                try:
                    print("There is no data for " + times[0].strftime('%m/%d/%Y %HZ') + " trying to retrieve data from the previous analysis at " + times[1].strftime('%m/%d/%Y %HZ'))
                    ds = xr.open_dataset(url_1, engine='netcdf4')
                    print("Data was successfully retrieved for " + times[1].strftime('%m/%d/%Y %HZ'))
                    strtime = times[1]
                    
                except Exception as b:
                        try:
                            print("There is no data for " + times[1].strftime('%m/%d/%Y %HZ') + " trying to retrieve data from the previous analysis at " + times[2].strftime('%m/%d/%Y %HZ'))
                            ds = xr.open_dataset(url_2, engine='netcdf4')
                            print("Data was successfully retrieved for " + times[2].strftime('%m/%d/%Y %HZ'))
                            strtime = times[2]
                            
                        except Exception as c:
                            try:
                                print("There is no data for " + times[2].strftime('%m/%d/%Y %HZ') + " trying to retrieve data from the previous analysis at " + times[3].strftime('%m/%d/%Y %HZ'))
                                ds = xr.open_dataset(url_3, engine='netcdf4')
                                print("Data was successfully retrieved for " + times[3].strftime('%m/%d/%Y %HZ'))
                                strtime = times[3]
                            except Exception as d:
    
                                try:
                                    print("There is no data for " + times[3].strftime('%m/%d/%Y %HZ') + " trying to retrieve data from the previous analysis at " + times[4].strftime('%m/%d/%Y %HZ'))
                                    ds = xr.open_dataset(url_4, engine='netcdf4')
                                    print("Data was successfully retrieved for " + times[4].strftime('%m/%d/%Y %HZ'))
                                    strtime = times[4]
                                    
                                except Exception as e:
                                    print("The latest dataset is over 4 hours old which isn't current. Please try again later.")
            
            temperature = ds['tmp2m']
            wind_speed = ds['wind10m']
            dewpoint = ds['dpt2m']
            temperature_k = units('kelvin') * temperature
            dewpoint_k = units('kelvin') * dewpoint
            relative_humidity = mpcalc.relative_humidity_from_dewpoint(temperature_k, dewpoint_k)
                
            lat = temperature['lat']
            lon = temperature['lon']
                
            lat_vals = lat[:].squeeze()
            lon_vals = lon[:].squeeze()

            # CONVERTS KELVIN TO FAHRENHEIT
            frac = 9/5
            temperature_f = (frac * (temperature - 273.15)) + 32
                

            # CONVERTS M/S TO MPH
            wind_speed = wind_speed * 2.23694
            
            ### This is what I want to ultimately plot with one shaded color where all these conditions hold true ###
            mask = (temperature_f.where(temperature_f >= 15) & relative_humidity.where(relative_humidity <= 25) & wind_speed.where(wind_speed >= 15))
                
            return lon_vals, lat_vals, strtime, mask

[DanielAdriaansen]

One approach I've used you can try:

# Mask is boolean
mask = (temperature_f.where(temperature_f >= 15) & \
        relative_humidity.where(relative_humidity <= 25) & \
        wind_speed.where(wind_speed >= 15))

# Create a DataArray full of zeros. This should have x-y dimensions and coordinates of the gridded data.
# I don't know what your coordinates and dims are, so I use x0 and y0 here.
plot_mask = xr.zeros_like(temperature,coords=temperature.coords['x0',y0'],dims=temperature.dims['x0','y0'])

# Switch the zeros to 1 where the mask is true (1)
plot_mask = xr.where(mask,1,plot_mask)

# Import matplotlib's colormap registry
from matplotlib.colors import ListedColormap

# Define a new colormap with two colors- white and the color for the hazard (here I use red).
# Using an alpha value of zero should make the white transparent so you can plot this on top of other fields,  I think.
mask_rgba = [(255/255,255/255,255/255,0),(255/255, 0/255, 0/255,1)]

# Register the colormap
mask_cmap = ListedColormap(mask_rgba)

# Plot the mask
ax.pcolormesh(lon,lat,plot_mask,transform=ccrs.PlateCarree(),cmap=mask_cmap,vmin=0.1)

[edrewitz]

Thanks for the reply! I tried the suggested method and got this error. I'm not sure what this error is supposed to mean.

[DanielAdriaansen]

Can you send the full traceback and the code that generated the error? In the sample code I sent, I assumed that:
temperature_f, relative_humidity, and wind_speed were all Xarray DataArray objects, but that may not be the case.

[edrewitz]

Thank you for the quick response! In the previous example, I was pulling data from NCEP NOMADS. In this case I tried doing the same thing with data from UCAR THREDDS since that data seems easier to work with in my experience. Here is the code:

def get_red_flag_warning_conditions(current_time):
        
            r"""
            THIS FUNCTION RETRIEVES THE LATEST 2.5KM X 2.5KM REAL TIME MESOSCALE ANALYSIS DATASETS FOR TEMPERATURE AND DEWPOINT
        
            THIS FUNCTION THEN CALCULATES A RELATIVE HUMIDITY DATASET USING METPY.CALC FROM THE TEMPERATURE AND DEWPOINT DATASETS
        
            IF THE DATASET FOR THE CURRENT TIME IS UNAVAILABLE THE FUNCTION WILL TRY TO RETURN THE MOST RECENT DATASET IN THE PAST 4 HOURS
        
            PYTHON PACKAGE DEPENDENCIES:
        
            1. SIPHON
            2. METPY
            3. DATETIME
        
            RETURNS:
        
            CURRENT RTMA DATASET FOR RELATIVE HUMIDITY
        
            COPYRIGHT (C) METEOROLOGIST ERIC J. DREWITZ 2024
        
            """
        
            times = []
        
            for i in range(0,5):
                new_time = current_time - timedelta(hours=i)
                times.append(new_time)
        
            try:
                rtma_cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/RTMA/CONUS_2p5km/RTMA_CONUS_2p5km_'+times[0].strftime('%Y%m%d_%H00')+'.grib2/catalog.xml')
                rtma_data = rtma_cat.datasets['RTMA_CONUS_2p5km_'+times[0].strftime('%Y%m%d_%H00')+'.grib2'].remote_access(use_xarray=True)
                rtma_data = rtma_data.metpy.parse_cf()
                rtma_temp = rtma_data['Temperature_Analysis_height_above_ground'].squeeze()
                rtma_dwpt = rtma_data['Dewpoint_temperature_Analysis_height_above_ground'].squeeze()
                rtma_wind = rtma_data['Wind_speed_Analysis_height_above_ground'].squeeze()
        
                rtma_rh = mpcalc.relative_humidity_from_dewpoint(rtma_temp, rtma_dwpt)
                print("Data retrieval for " + times[0].strftime('%m/%d/%Y %H00 UTC') + " is successful")

                time = times[0]
                
                return rtma_rh *100, rtma_wind * 2.23694, time
                
            except Exception as e:
        
                print("Relative Humidity Data is unavailiable for "+times[0].strftime('%m/%d/%Y %H00 UTC')+ "\nWill try to download the most recent dataset from "+times[1].strftime('%m/%d/%Y %H00 UTC'))
                
                try:
                    rtma_cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/RTMA/CONUS_2p5km/RTMA_CONUS_2p5km_'+times[1].strftime('%Y%m%d_%H00')+'.grib2/catalog.xml')
                    rtma_data = rtma_cat.datasets['RTMA_CONUS_2p5km_'+times[1].strftime('%Y%m%d_%H00')+'.grib2'].remote_access(use_xarray=True)
                    rtma_data = rtma_data.metpy.parse_cf()
                    rtma_temp = rtma_data['Temperature_Analysis_height_above_ground'].squeeze()
                    rtma_dwpt = rtma_data['Dewpoint_temperature_Analysis_height_above_ground'].squeeze()
                    rtma_wind = rtma_data['Wind_speed_Analysis_height_above_ground'].squeeze()
            
                    rtma_rh = mpcalc.relative_humidity_from_dewpoint(rtma_temp, rtma_dwpt)
                    print("Data retrieval for " + times[1].strftime('%m/%d/%Y %H00 UTC') + " is successful")

                    time = times[1]
                    
                    return rtma_rh *100, rtma_wind, time
          
                except Exception as a:
        
                    print("Relative Humidity data is unavailiable for "+times[1].strftime('%m/%d/%Y %H00 UTC')+ "\nWill try to download the most recent dataset from "+times[2].strftime('%m/%d/%Y %H00 UTC'))
                    
                    try:
                        rtma_cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/RTMA/CONUS_2p5km/RTMA_CONUS_2p5km_'+times[2].strftime('%Y%m%d_%H00')+'.grib2/catalog.xml')
                        rtma_data = rtma_cat.datasets['RTMA_CONUS_2p5km_'+times[2].strftime('%Y%m%d_%H00')+'.grib2'].remote_access(use_xarray=True)
                        rtma_data = rtma_data.metpy.parse_cf()
                        rtma_temp = rtma_data['Temperature_Analysis_height_above_ground'].squeeze()
                        rtma_dwpt = rtma_data['Dewpoint_temperature_Analysis_height_above_ground'].squeeze()
                        rtma_wind = rtma_data['Wind_speed_Analysis_height_above_ground'].squeeze()
                
                        rtma_rh = mpcalc.relative_humidity_from_dewpoint(rtma_temp, rtma_dwpt)
                        print("Data retrieval for " + times[2].strftime('%m/%d/%Y %H00 UTC') + " is successful")

                        time = times[2]
                        
                        return rtma_rh *100, rtma_wind, time
        
                    except Exception as b:
                                    
                        print("Relative Humidity data is unavailiable for "+times[2].strftime('%m/%d/%Y %H00 UTC')+ "\nWill try to download the most recent dataset from "+times[3].strftime('%m/%d/%Y %H00 UTC'))
                        
                        try:
                            rtma_cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/RTMA/CONUS_2p5km/RTMA_CONUS_2p5km_'+times[3].strftime('%Y%m%d_%H00')+'.grib2/catalog.xml')
                            rtma_data = rtma_cat.datasets['RTMA_CONUS_2p5km_'+times[3].strftime('%Y%m%d_%H00')+'.grib2'].remote_access(use_xarray=True)
                            rtma_data = rtma_data.metpy.parse_cf()
                            rtma_temp = rtma_data['Temperature_Analysis_height_above_ground'].squeeze()
                            rtma_dwpt = rtma_data['Dewpoint_temperature_Analysis_height_above_ground'].squeeze()
                            rtma_wind = rtma_data['Wind_speed_Analysis_height_above_ground'].squeeze()
                    
                            rtma_rh = mpcalc.relative_humidity_from_dewpoint(rtma_temp, rtma_dwpt)
                            print("Data retrieval for " + times[3].strftime('%m/%d/%Y %H00 UTC') + " is successful")

                            time = times[3]
                            
                            return rtma_rh *100, rtma_wind, time
        
                        except Exception as c:
                                    
                            print("Relative Humidity data is unavailiable for "+times[3].strftime('%m/%d/%Y %H00 UTC')+ "\nWill try to download the most recent dataset from "+times[4].strftime('%m/%d/%Y %H00 UTC'))
                            
                            try:
                                rtma_cat = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/RTMA/CONUS_2p5km/RTMA_CONUS_2p5km_'+times[4].strftime('%Y%m%d_%H00')+'.grib2/catalog.xml')
                                rtma_data = rtma_cat.datasets['RTMA_CONUS_2p5km_'+times[4].strftime('%Y%m%d_%H00')+'.grib2'].remote_access(use_xarray=True)
                                rtma_data = rtma_data.metpy.parse_cf()
                                rtma_temp = rtma_data['Temperature_Analysis_height_above_ground'].squeeze()
                                rtma_dwpt = rtma_data['Dewpoint_temperature_Analysis_height_above_ground'].squeeze()
                                rtma_wind = rtma_data['Wind_speed_Analysis_height_above_ground'].squeeze()
                        
                                rtma_rh = mpcalc.relative_humidity_from_dewpoint(rtma_temp, rtma_dwpt)
                                print("Data retrieval for " + times[4].strftime('%m/%d/%Y %H00 UTC') + " is successful")

                                time = times[4]
                                
                                return rtma_rh *100, rtma_wind, time
                                       
                        
                            except Exception as k:
                                print("WARNING: Latest dataset is more than 4 hours old.\nQuitting - Please try again later.")
        
                                return None

Then in a separate jupyter notebook here is what I did(the function I call is the function shown above):

import data_access as da
from NWS_Generic_Forecast_Graphics import standard

local, utc = standard.plot_creation_time()
mask = 30000

rtma_rh, rtma_wind, time = da.UCAR_THREDDS_SERVER_OPENDAP_Downloads.CONUS.get_red_flag_warning_conditions(utc)

mask = (rtma_rh.where(rtma_rh <= 25) & rtma_wind.where(rtma_wind >= 15))

Here is the full traceback:

---

TypeError Traceback (most recent call last)
Cell In[7], line 1
----> 1 mask = (rtma_rh.where(rtma_rh <= 25) & rtma_wind.where(rtma_wind >= 15))

File ~\miniconda3\lib\site-packages\xarray\core_typed_ops.py:266, in DataArrayOpsMixin.and(self, other)
265 def and(self, other: DaCompatible) -> Self:
--> 266 return self.binary_op(other, operator.and)

File ~\miniconda3\lib\site-packages\xarray\core\dataarray.py:4691, in DataArray._binary_op(self, other, f, reflexive)
4687 other_variable_or_arraylike: DaCompatible = getattr(other, "variable", other)
4688 other_coords = getattr(other, "coords", None)
4690 variable = (
-> 4691 f(self.variable, other_variable_or_arraylike)
4692 if not reflexive
4693 else f(other_variable_or_arraylike, self.variable)
4694 )
4695 coords, indexes = self.coords._merge_raw(other_coords, reflexive)
4696 name = self._result_name(other)

File ~\miniconda3\lib\site-packages\xarray\core_typed_ops.py:542, in VariableOpsMixin.and(self, other)
541 def and(self, other: VarCompatible) -> Self | T_DataArray:
--> 542 return self.binary_op(other, operator.and)

File ~\miniconda3\lib\site-packages\xarray\core\variable.py:2411, in Variable._binary_op(self, other, f, reflexive)
2408 attrs = self._attrs if keep_attrs else None
2409 with np.errstate(all="ignore"):
2410 new_data = (
-> 2411 f(self_data, other_data) if not reflexive else f(other_data, self_data)
2412 )
2413 result = Variable(dims, new_data, attrs=attrs)
2414 return result

TypeError: operand type(s) all returned NotImplemented from array_ufunc(<ufunc 'bitwise_and'>, 'call', <Quantity([[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]
...
[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]], 'dimensionless')>, array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]], dtype=float32)): 'Quantity', 'ndarray'

I also checked by looking at the individual parameter data arrays after downloading the data before trying to apply the mask to verify that there is actual data rather than [nan, nan, nan] and there was so that adds to my uncertainty of this error. I also got this same error in the prior example with the data from NCEP NOMADS. Thank you so much for your help, I appreciate it!

[DanielAdriaansen]

Try this:
mask = (rtma_rh <= 25) & (rtma_wind >= 15)

I think it's my oops for suggesting where() when writing the pseudo-code!

[edrewitz]

Thank you so much for your help! My plot works now. I will post the rest of my code here to show how I finished this off in case anyone else has this same question in the future so they have an example to start from. I changed the mask to reflect lower-48 red flag warning criteria in this example. When running this originally I got an error since with the current weather pattern there were no red flag conditions anywhere in the area. I then solved the issue with a simple try and except block with the ax.pcolormesh line in the try block and a pass in the except block to allow for a blank map when none of these conditions are met anywhere on the map.

Here is the example code:

        def plot_red_flag_criteria_based_on_sustained_winds_test(western_bound, eastern_bound, southern_bound, northern_bound, central_longitude, central_latitude, first_standard_parallel, second_standard_parallel, fig_x_length, fig_y_length, color_table_shrink, signature_x_position, signature_y_position, title_font_size, signature_font_size, colorbar_label_font_size, colorbar_pad):

            r'''
            THIS FUNCTION CREATES A CUSTOMIZED PLOT OF THE 2.5KM X 2.5KM REAL TIME MESOSCALE ANALYSIS DATA FOR ANY AREA INSIDE OF CONUS. THIS PLOT COMPARES THE AREAS OF RED-FLAG RELATIVE HUMIDITY CRITERIA WITH RED-FLAG WIND CRITERIA BASED ON SUSTAINED WINDS. 

            (C) METEOROLOGIST ERIC J. DREWITZ 2024
            
            '''
            colorbar_label_font_size = colorbar_label_font_size

            colorbar_pad = colorbar_pad

            local_time, utc_time = standard.plot_creation_time()

            rtma_rh, rtma_wind, rtma_time = da.UCAR_THREDDS_SERVER_OPENDAP_Downloads.CONUS.get_red_flag_warning_conditions(utc_time)

            rtma_wind = rtma_wind * 2.23694
            
            mapcrs = ccrs.LambertConformal(central_longitude=central_longitude, central_latitude=central_latitude, standard_parallels=(first_standard_parallel,second_standard_parallel))
            datacrs = ccrs.PlateCarree()

            mask = (rtma_rh <= 15) & (rtma_wind >= 25)
            lon = mask['longitude']
            lat = mask['latitude']

            plot_proj = mask.metpy.cartopy_crs

            fig = plt.figure(figsize=(fig_x_length, fig_y_length))

            ax = fig.add_subplot(1, 1, 1, projection=plot_proj)
            ax.set_extent((western_bound, eastern_bound, southern_bound, northern_bound), crs=ccrs.PlateCarree())
            ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
            ax.add_feature(cfeature.STATES, linewidth=2, edgecolor='blue', zorder=3)
            ax.add_feature(USCOUNTIES, linewidth=1.5, zorder=2)
            
            
            # Import matplotlib's colormap registry
            from matplotlib.colors import ListedColormap
            
            mask_rgba = [(255/255,255/255,255/255,0),(255/255, 0/255, 0/255,1)]
            
            # Register the colormap
            mask_cmap = ListedColormap(mask_rgba)
            
            # Plot the mask
            try:
                ax.pcolormesh(lon,lat,mask,transform=ccrs.PlateCarree(),cmap=mask_cmap,vmin=0.1)

            except Exception as e:
                pass


            plt.title("Red-Flag Warning Conditions (RH <= 15% and Wind Speed >= 25 MPH)\nValid: " + rtma_time.strftime('%m/%d/%Y %HZ') + "\nImage Created: " + utc_time.strftime('%m/%d/%Y %H:%MZ'), fontsize=title_font_size, fontweight='bold')
            
            ax.text(signature_x_position, signature_y_position, "Plot Created With FireWxPy (C) Eric J. Drewitz 2024\nData Source: thredds.ucar.edu", fontsize=signature_font_size, fontweight='bold', horizontalalignment='center',
           verticalalignment='bottom', transform=ax.transAxes)

            return fig