Alexander :
Attached is a little "command line" python script I wrote. It uses
xarray and metpy.declarative.
You can make a MSL chart by running it as :
python xarray-declarative-2.py --nwpfield Pressure_reduced_to_MSL_msl -a opcahsf
... or "500 mb" chart with :
python xarray-declarative-2.py -f Geopotential_height_isobaric -l '500
hPa' -a opcahsf
The default is to make a chart closest to the current time. You can
specify the "--reftime" with a command line option. If you specify in
the future, it will use a forecast. For example :
python xarray-declarative-2.py --nwpfield Pressure_reduced_to_MSL_msl
-a opcahsf --reftime 2021-08-23T00:00
It uses the thredds "Best" time to get the most "best" reftime and
valid time. There is a poster that explains this ... which I can't
seem to find right now.
Sorry ... I don't print the reftime / valid time on the plot. I can
try to figure that out if you want.
Thanks,
Ken
On Tue, Aug 24, 2021 at 1:00 PM Alexander Martínez
<almmartinezme@xxxxxxxxx> wrote:
>
> Hello.I made an example of the surface chart linking
> https://github.com/Unidata/python-training/blob/master/pages/gallery/HILO_Symbol_Plot.ipynb,
> but the Thredds catalog has data from 2014 and not data from August 2021. I
> also have problems with the pressure levels because the colors of the high
> and low pressures are not displayed.
>
> Could you please show me an example of a surface chart for the area, but with
> data from August 23?
>
> Thank you very much for your collaboration and I remain attentive to your
> instructions.
>
>
>
> --
> Alexander M. Martínez Mercado
> Lic. en Química
> Universidad Distrital Francisco José de Caldas
> MSc en Ciencias - Meteorología
> Universidad Nacional de Colombia
> Docente MT UECCI
> OSPA - IDEAM
> _______________________________________________
> NOTE: All exchanges posted to Unidata maintained email lists are
> recorded in the Unidata inquiry tracking system and made publicly
> available through the web. Users who post to any of the lists we
> maintain are reminded to remove any personal information that they
> do not want to be made public.
>
>
> python-users mailing list
> python-users@xxxxxxxxxxxxxxxx
> For list information, to unsubscribe, or change your membership options,
> visit: https://www.unidata.ucar.edu/mailing_lists/
#!/usr/bin/python3
# plot with "metpy declarative"
# usage :
# python xarray-declarative-2.py -f Pressure_reduced_to_MSL_msl -a opcahsf
# python xarray-declarative-2.py -f Geopotential_height_isobaric -l '500 hPa'
# on Fedora :
# dnf install python3-xarray python3-cartopy
# pip install metpy
#############################################################
import xarray
import metpy
def proc0(args1):
#print('=== args1 : ', args1)
args0 = vars(args1) # usings vars() to create dict
print('=== args0 : ', args0)
########################
if args0['level'] is None:
lev1 = None
else:
s0 = args0['level'].split()
lev1 = int(s0[0]) * metpy.units.units(s0[1])
########################
import datetime
if args0['reftime'] is not None:
reftime0 = datetime.datetime.strptime(args0['reftime'],
'%Y-%m-%dT%H:%M')
else:
reftime0 = datetime.datetime.utcnow()
vtime0 = reftime0 + datetime.timedelta(hours=int(args0['valid'])) #
forecast time / valid time
########################
n1 = args0['nwpmodel']
url0 = f'https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/{n1}/Best'
xdsnwp0 = xarray.open_dataset(url0)
#print(xdsnwp0)
from metpy.plots import declarative
panel0 = declarative.MapPanel()
######################################### plot NWP data
# Set Contour Plot Parameters
#contour0 = declarative.FilledContourPlot()
contour0 = declarative.ContourPlot()
contour0.data = xdsnwp0
contour0.time = vtime0
contour0.field = args0['nwpfield']
contour0.level = lev1
contour0.clabels = True
contour0.linecolor = 'blue'
contour0.linewidth = 1
v0 = contour0.data[contour0.field]
print(v0)
#if args0['unit'] is not None:
#v0.metpy.convert_units(args0['unit']) # XXX broken ?? slow
panel0.title += f'{n1} '
panel0.title += f'{v0.name}'
if contour0.level is not None:
panel0.title += f'@{contour0.level}'
panel0.title += " [{}]".format(v0.attrs["units"])
#panel0.title += '\n' ###################
#panel0.title += xarray.core.formatting.format_item(v0.reftime.values)
#########################################
panel0.area = args0['area'] # XXX default : max of all plots ?? see
metpy/plots/declarative.py : _areas
if args0['proj'] is not None:
panel0.projection = args0['proj'] # XXX default : proj of first plot
import cartopy.feature
panel0.layers = [cartopy.feature.COASTLINE]#(edgecolor='white')]
panel0.plots = [contour0]
pc0 = declarative.PanelContainer()
pc0.panels = [panel0]
pc0.show()
#pc0.save(f'fn_{dt:%Y%m%d_%H%MZ}.{ext}', dpi=200, bbox_inches='tight')
# del pc0
def main():
import argparse
parser = argparse.ArgumentParser(description="get point data")
parser.add_argument("-m", "--nwpmodel", default = 'GFS/Global_onedeg',
help="NWP model")
parser.add_argument("-f", "--nwpfield", default =
'Pressure_reduced_to_MSL_msl', help="field")
parser.add_argument("-a", "--area", default = None, type = str, help="area")
parser.add_argument("-p", "--proj", help="projection", default = None, type
= str)
#parser.add_argument("-u", "--unit", default = None , help="plot units",
type = str)
# -u knot
# -u hPa
parser.add_argument("-l", "--level", default = None , help="level", type =
str)
# -l '500 hPa'
parser.add_argument("-t", "--reftime", default = None, help="ref time :
%Y-%m-%dT%H:%M")
parser.add_argument("-v", "--valid", default = 0, type = int, help="valid
time (hours)")
proc0(parser.parse_args())
import sys
if __name__ == '__main__':
sys.exit(main())
