Surface map of METAR data

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

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

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

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

# 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):
    time_1 = datetime.now()
    offset = timedelta(hours = 1) # Worldwide METARS typically don't come in until ~15 minutes past the hour
    time_1 = time_1 - offset
else:
    year = 2026
    month = 6
    day = 18
    hour = 18
    time_1 = datetime(year, month, day, hour)

time_0 = time_1 - timedelta(minutes=15)
YYYY_0 = time_0.strftime("%Y")
YYYY_1 = time_1.strftime("%Y")
date_string = time_1.strftime("%Y-%m-%d %H00 UTC")
print(time_0, time_1)
# Handle edge case when the two hours straddle the end/beginning of a yearw
if (YYYY_0 == YYYY_1):
    URLs = [f'https://data.source.coop/dynamical/asos-parquet/year={YYYY_1}/data.parquet']
else:
    URLs = [f'https://data.source.coop/dynamical/asos-parquet/year={YYYY_0}/data.parquet',
            f'https://data.source.coop/dynamical/asos-parquet/year={YYYY_1}/data.parquet']

df = duckdb.execute("""
    SELECT *
    FROM read_parquet($1, hive_partitioning=true)
    WHERE 
---      country = 'FR' AND
---    station = 'ENGM' AND
    valid BETWEEN $2 AND $3
    ORDER BY country
""", [URLs, time_0, time_1]).fetchdf()

extend = 1
df2 = df.query('latitude >= (@latS - @extend) & latitude <= (@latN + @extend) & longitude  >= (@lonW - @extend) & longitude <= (@lonE + @extend)')

df2

df2.columns

columnSubset = ['station', 'valid', 'latitude', 'longitude', 'elevation', 'tmpc', 'dwpc', 'mslp',
                'sknt', 'drct','alti','vsby','p01i']

#columnSubset = ['station', 'valid', 'latitude', 'longitude', 'elevation', 'tmpc', 'dwpc', 'mslp',
#                'sknt', 'drct','alti','vsby','p01i','VSBY','CHC1', 'CHC2', 'CHC3','CTYH', 'CTYM', 'CTYL']

df3 = df2[columnSubset]

df3

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

df4

data = df4

lats = data['latitude'].values
lons = data['longitude'].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['mslp'].values * units('hPa')

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

stid = data['station']

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

proj_data = ccrs.PlateCarree()
spacing = 150000 # meters; change depending on the map area
xy = proj_map.transform_points(proj_data, lons, lats)
mask = reduce_point_density(xy, spacing)

# 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=(20,12), dpi=dpi)

ax = fig.add_subplot(1, 1, 1, projection=proj_map) # 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())

# Slightly expand 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
w_ext = -0.5
e_ext = 0.5
s_ext = -1
n_ext = 1
ax.set_extent ((lonW + w_ext,lonE + e_ext,latS + s_ext, latN + n_ext), crs=proj_data)

#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)

# Remove any text that appears outside the borders of the Axes

for obj in ax.findobj(mtext.Text):
    obj.set_clip_on(True)
plotTitle = (f"Sfc Map valid at:  {date_string}")
ax.set_title (plotTitle);