[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[python #OTR-911008]: problem with Example MSLP and 1000-500 hPa Thickness with High and Low Symbols



Hi! Glad I was able to point you in the right direction, and thanks for 
checking back in.

The issue here is that that the NARR record was no longer ingested by NCEI as 
of 2014 (see: 
https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/north-american-regional-reanalysis-narr).
 If you need newer NARR records, you can check out NCAR's research data archive 
(RDA) record for the NARR, available at https://rda.ucar.edu/datasets/ds608.0/ 
or over at NOAA PSL at https://psl.noaa.gov/data/narr/ where you can find 
access through FTP or their own THREDDS data server (TDS).

You can also explore using this code with other products, data, and output! I 
hope this helps, and please let me know if I can help out with the next steps.


All the best,

Drew


> Drew,
> 
> thank you very much.
> I switched to the new ntebook code you provided.
> But now I get an ugly "Segmentation fault".
> 
> Is there a debug switch I can switch on to get more info?
> 
> I'm not very skilled and don't want to bother you with my problems.
> The idea is to find some code to plot the recent MSLP and high and low
> spots on top of a satpy generated background composite form MSG4 sat data.
> If I search for some code I always ended up to MetPy.
> 
> Maybe you have some idea howto provide this goal I want to achieve?
> I tried some actual data for
> 
> dattim = datetime(2020, 6, 11, 0)
> 
> (segfault too but gives:
> 
> requests.exceptions.HTTPError: Error accessing
> https://www.ncei.noaa.gov/thredds/ncss/model-narr-a-files/202006/20200611/narr-a_221_20200611_0000_000.grb/dataset.xml
> Server Error (404:
> model-narr-a-files/202006/20200611/narr-a_221_20200611_0000_000.grb)
> 
> but more recent datait seems not available at this source.
> 
> What is the idea: to plot most recent data on top of the recent sat pic).
> Is this code/souerce for MSLP then the right place to achive my idea!?
> 
> 
> Thanks in advance and for you patience!
> 
> 
> 
> This is my code wich segfaults (it's a script (not inside notbook) and I
> would save the output into a file):
> 
> from datetime import datetime
> 
> import cartopy.crs as ccrs
> import cartopy.feature as cfeature
> import matplotlib.pyplot as plt
> from metpy.units import units
> from netCDF4 import num2date
> import numpy as np
> from scipy.ndimage import gaussian_filter
> from siphon.ncss import NCSS
> 
> def plot_maxmin_points(lon, lat, data, extrema, nsize, symbol, color='k',
> plotValue=True, transform=None):
> """
> This function will find and plot relative maximum and minimum for a
> 2D grid. The function
> can be used to plot an H for maximum values (e.g., High pressure)
> and an L for minimum
> values (e.g., low pressue). It is best to used filetered data to
> obtain  a synoptic scale
> max/min value. The symbol text can be set to a string value and
> optionally the color of the
> symbol and any plotted value can be set with the parameter color
> lon = plotting longitude values (2D)
> lat = plotting latitude values (2D)
> data = 2D data that you wish to plot the max/min symbol placement
> extrema = Either a value of max for Maximum Values or min for
> Minimum Values
> nsize = Size of the grid box to filter the max and min values to
> plot a reasonable number
> symbol = String to be placed at location of max/min value
> color = String matplotlib colorname to plot the symbol (and
> numerica value, if plotted)
> plot_value = Boolean (True/False) of whether to plot the numeric
> value of max/min point
> The max/min symbol will be plotted on the current axes within the
> bounding frame
> (e.g., clip_on=True)
> """
> from scipy.ndimage.filters import maximum_filter, minimum_filter
> 
> if (extrema == 'max'):
> data_ext = maximum_filter(data, nsize, mode='nearest')
> elif (extrema == 'min'):
> data_ext = minimum_filter(data, nsize, mode='nearest')
> else:
> raise ValueError('Value for hilo must be either max or min')
> 
> mxy, mxx = np.where(data_ext == data)
> 
> for i in range(len(mxy)):
> ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]], symbol,
> color=color, size=24,
> clip_on=True, horizontalalignment='center',
> verticalalignment='center',
> transform=transform)
> ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]],
> '\n' + str(np.int(data[mxy[i], mxx[i]])),
> color=color, size=12, clip_on=True, fontweight='bold',
> horizontalalignment='center', verticalalignment='top',
> transform=transform)
> 
> dattim = datetime(1999, 1, 3, 0)
> #dattim = datetime(2020, 6, 11, 0)
> ncss =
> NCSS('https://www.ncei.noaa.gov/thredds/ncss/model-narr-a-files/{0:%Y%m}/{0:%Y%m%d}/'
> 'narr-a_221_{0:%Y%m%d}_{0:%H}00_000.grb'.format(dattim))
> 
> query = ncss.query()
> query.all_times().variables('Pressure_reduced_to_MSL_msl',
> 
> 'Geopotential_height_isobaric').add_lonlat().accept('netcdf')
> data = ncss.get_data(query)
> 
> # Grab pressure levels
> plev = list(data.variables['isobaric1'][:])
> 
> # Grab lat/lons and make all lons 0-360
> lats = data.variables['lat'][:]
> lons = data.variables['lon'][:]
> lons[lons < 0] = 360 + lons[lons < 0]
> 
> # Grab valid time and get into datetime format
> time = data['time2']
> vtime = num2date(time[:], units=time.units)
> 
> # Grab MSLP and smooth, use MetPy Units module for conversion
> emsl_var = data.variables['Pressure_reduced_to_MSL_msl']
> EMSL = units.Quantity(emsl_var[:], emsl_var.units).to('hPa')
> #EMSL = units.Quantity(data.variables['Pressure_reduced_to_MSL_msl'][:],
> units.Pa).to('hPa')
> mslp = gaussian_filter(EMSL[0], sigma=3.0)
> 
> # Grab pressure level data
> hght_1000 = data.variables['Geopotential_height_isobaric'][0,
> plev.index(1000)]
> hght_500 = data.variables['Geopotential_height_isobaric'][0,
> plev.index(500)]
> 
> # Calculate and smooth 1000-500 hPa thickness
> thickness_1000_500 = gaussian_filter(hght_500 - hght_1000, sigma=3.0)
> 
> # Set projection of map display
> mapproj = ccrs.LambertConformal(central_latitude=45.,
> central_longitude=-100.)
> 
> # Set projection of data
> dataproj = ccrs.PlateCarree()
> 
> # Grab data for plotting state boundaries
> states_provinces = cfeature.NaturalEarthFeature(
> category='cultural',
> name='admin_1_states_provinces_lakes',
> scale='50m',
> facecolor='none')
> 
> fig = plt.figure(1, figsize=(17., 11.))
> ax = plt.subplot(111, projection=mapproj)
> 
> # Set extent and plot map lines
> ax.set_extent([-145., -70, 20., 60.], ccrs.PlateCarree())
> ax.coastlines('50m', edgecolor='black', linewidth=0.75)
> ax.add_feature(states_provinces, edgecolor='black', linewidth=0.5)
> 
> # Plot thickness with multiple colors
> clevs = (np.arange(0, 5400, 60),
> np.array([5400]),
> np.arange(5460, 7000, 60))
> colors = ('tab:blue', 'b', 'tab:red')
> kw_clabels = {'fontsize': 11, 'inline': True, 'inline_spacing': 5,
> 'fmt': '%i',
> 'rightside_up': True, 'use_clabeltext': True}
> for clevthick, color in zip(clevs, colors):
> cs = ax.contour(lons, lats, thickness_1000_500, levels=clevthick,
> colors=color,
> linewidths=1.0, linestyles='dashed',
> transform=dataproj)
> plt.clabel(cs, **kw_clabels)
> 
> # Plot MSLP
> clevmslp = np.arange(800., 1120., 4)
> cs2 = ax.contour(lons, lats, mslp, clevmslp, colors='k', linewidths=1.25,
> linestyles='solid', transform=dataproj)
> plt.clabel(cs2, **kw_clabels)
> 
> # Use definition to plot H/L symbols
> plot_maxmin_points(lons, lats, mslp, 'max', 50, symbol='H', color='b',
> transform=dataproj)
> plot_maxmin_points(lons, lats, mslp, 'min', 25, symbol='L', color='r',
> transform=dataproj)
> 
> # Put on some titles
> plt.title('MSLP (hPa) with Highs and Lows, 1000-500 hPa Thickness (m)',
> loc='left')
> plt.title('VALID: {}'.format(vtime[0]), loc='right')
> 
> #plt.save('/tmp/test.png')
> 
> 

Ticket Details
===================
Ticket ID: OTR-911008
Department: Support Python
Priority: Low
Status: Closed
===================
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata 
inquiry tracking system and then made publicly available through the web.  If 
you do not want to have your interactions made available in this way, you must 
let us know in each email you send to us.