Surface map of METAR data

from datetime import datetime,timedelta
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import cartopy.crs as ccrs
import cartopy.feature as cfeature

from metpy.calc import wind_components, reduce_point_density
from metpy.units import units
from metpy.plots import StationPlot
from metpy.plots.wx_symbols import current_weather, sky_cover, wx_code_map

# Load in a collection of functions that process GEMPAK weather conditions and cloud cover data.
%run /kt11/ktyle/python/metargem.py

# Use the current time, or set your own for a past time.
# Set current to False if you want to specify a past time.

nowTime = datetime.now()

current = True
#current = False
if (current):
    validTime = datetime.now()
    year = validTime.year
    month = validTime.month
    day = validTime.day
    hour = validTime.hour
    validTime = datetime(year, month, day, hour)
    offset = timedelta(hours = 1)
    validTime = validTime - offset
else:
    year = 2010
    month = 2
    day = 5
    hour = 23
    validTime = datetime(year, month, day, hour)

timeStr = validTime.strftime("%Y-%m-%d %H UTC")
timeStr2 = validTime.strftime("%Y%m%d%H")
YYMMDDHH = validTime.strftime("%y%m%d%H")
print(timeStr)
print(validTime)

metarFile = f'/ktyle_rit/scripts/sflist2/complete/{YYMMDDHH}.csv'
df = pd.read_csv(metarFile, sep='\\s+')

df

latN = 55
latS = 20
lonW = -125
lonE = -60

cLat = (latN + latS)/2
cLon = (lonW + lonE )/2

extend = .01
df2 = df.query('SLAT >= (@latS - @extend) & SLAT <= (@latN + @extend) & SLON  >= (@lonW - @extend) & SLON <= (@lonE + @extend)')

df2

columnSubset = ['STN', 'YYMMDD/HHMM', 'SLAT', 'SLON', 'SELV', 'TMPC', 'DWPC', 'PMSL',
                'SKNT', 'DRCT','ALTI','WNUM','VSBY','CHC1', 'CHC2', 'CHC3','CTYH', 'CTYM', 'CTYL']

df3 = df2[columnSubset]

df3

df4 = df3.replace(-9999.0, np.nan)

df4

data = df4

data.loc[data['WNUM'] =='********', ['WNUM']] = '-9999.00'
data['WNUM'] = data['WNUM'].astype('float16')

lats = data['SLAT'].values
lons = data['SLON'].values
tair = (data['TMPC'].values * units ('degC')).to('degF')
dewp = (data['DWPC'].values * units ('degC')).to('degF')
altm = (data['ALTI'].values * units('inHg')).to('mbar')
slp = data['PMSL'].values * units('hPa')

# Convert wind to components
u, v = wind_components(data['SKNT'].values * units.knots, data['DRCT'].values * units.degree)

# replace missing wx codes or those >= 100 with 0, and convert to MetPy's present weather code
wnum = (np.nan_to_num(data['WNUM'].values,True).astype(int))
convert_wnum (wnum)

# Need to handle missing (NaN) and convert to proper code
chc1 = (np.nan_to_num(data['CHC1'].values,True).astype(int))
chc2 = (np.nan_to_num(data['CHC2'].values,True).astype(int))
chc3 = (np.nan_to_num(data['CHC3'].values,True).astype(int))
cloud_cover = calc_clouds(chc1, chc2, chc3)

# For some reason station id's come back as bytes instead of strings. This line converts the array of station id's into an array of strings.
#stid = np.array([s.decode() for s in data['station']])
stid = data['STN']

# Uncomment whichever projection you want (or add your own!)
#proj = ccrs.Stereographic(central_longitude=cLon, central_latitude=cLat)
proj = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)

# Project points so that we're filtering based on the way the stations are laid out on the map

xy = proj.transform_points(ccrs.PlateCarree(), lons, lats)

# Reduce point density so that there's only one point within a circle whose distance is specified in meters.
# This value will need to change depending on how large of an area you are plotting.
density = 150000
mask = reduce_point_density(xy, density)

# Set up a plot with map features
# First set dpi ("dots per inch") - higher values will give us a less pixelated final figure.
dpi = 125

fig = plt.figure(figsize=(24,18), dpi=dpi)

ax = fig.add_subplot(1, 1, 1, projection=proj) # We use the projection defined above
ax.set_facecolor(cfeature.COLORS['water'])

land_mask = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])
lake_mask = cfeature.NaturalEarthFeature('physical', 'lakes', '50m',
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['water'])
state_borders = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lakes',
                                         scale='50m', facecolor='none')

ax.add_feature(land_mask)
ax.add_feature(lake_mask)
ax.add_feature(state_borders, linestyle='solid', edgecolor='black')

# Slightly reduce the extent of the map as compared to the subsetted data region; this helps eliminate data from being plotted beyond the frame of the map
ax.set_extent ((lonW+2,lonE-4,latS+1,latN-2), crs=ccrs.PlateCarree())
ax.set_extent ((lonW+0.2,lonE-0.2,latS+.1,latN-.1), crs=ccrs.PlateCarree())
ax.set_extent ((lonW,lonE,latS,latN), crs=ccrs.PlateCarree())

#If we wanted to add grid lines to our plot:
#ax.gridlines()

# Create a station plot pointing to an Axes to draw on as well as the location of points
stationplot = StationPlot(ax, lons[mask], lats[mask], transform=ccrs.PlateCarree(),
                          fontsize=8)

stationplot.plot_parameter('NW', tair[mask], color='red', fontsize=10)
stationplot.plot_parameter('SW', dewp[mask], color='darkgreen', fontsize=10)

# Below, we are using a custom formatter to control how the sea-level pressure
# values are plotted. This uses the standard trailing 3-digits of the pressure value
# in tenths of millibars.
stationplot.plot_parameter('NE', slp[mask], color='purple', formatter=lambda v: format(10 * v, '.0f')[-3:])



stationplot.plot_symbol('C', cloud_cover[mask], sky_cover)
stationplot.plot_symbol('W', wnum[mask], current_weather,color='blue',fontsize=12)
stationplot.plot_text((2, 0),stid[mask], color='gray')
#zorder - Higher value zorder will plot the variable on top of lower value zorder. This is necessary for wind barbs to appear. Default is 1.
stationplot.plot_barb(u[mask], v[mask],zorder=2)

plotTitle = (f"Sfc Map valid at:  {timeStr}")
ax.set_title (plotTitle);

# In order to see the entire figure, type the name of the figure object below.
fig

figName = (f'{timeStr2}_sfmap.png')
fig.savefig(figName)