Cartopy 2: NYS Mesonet Data

In this notebook, we’ll use Cartopy, Matplotlib, Datetime, and Pandas to visualize data from the New York State Mesonet, headquartered right here at UAlbany.

import matplotlib.pyplot as plt
import pandas as pd
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from datetime import datetime

Create a regional map, centered over NYS, and add in some geographic features.

Be patient: this may take a minute or so to plot, depending on the resolution of the Natural Earth shapefile features you are adding!

Create the figure

For a quick look, let’s just choose the coarsest (110,000,000:1) Natural Earth shapefiles set.

# Set the domain for defining the plot region.
latN = 45.2
latS = 40.2
lonW = -80.0
lonE = -71.5
cLat = (latN + latS)/2
cLon = (lonW + lonE )/2


proj = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)

res = '110m' # Coarsest and quickest to display; other options are '10m' (slowest) and '50m'.

fig = plt.figure(figsize=(11,8.5),dpi=125)
ax = plt.subplot(1,1,1,projection=proj)
ax.set_extent ([lonW,lonE,latS,latN])
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res));

Plot some data on the map. We’ll use Pandas to read in the file containing the most recent NYS Mesonet obs.

df = pd.read_csv('http://www.atmos.albany.edu/products/nysm/nysm_latest.csv')

View the first and last five lines of this DataFrame

df

	station	time	temp_2m [degC]	temp_9m [degC]	relative_humidity [percent]	precip_incremental [mm]	precip_local [mm]	precip_max_intensity [mm/min]	avg_wind_speed_prop [m/s]	max_wind_speed_prop [m/s]	...	soil_temp_05cm [degC]	soil_temp_25cm [degC]	soil_temp_50cm [degC]	soil_moisture_05cm [m^3/m^3]	soil_moisture_25cm [m^3/m^3]	soil_moisture_50cm [m^3/m^3]	lat	lon	elevation	name
0	ADDI	2025-03-27 15:30:00	2.2	1.9	41.8	0.0	0.00	0.0	3.5	5.5	...	1.0	2.6	3.5	0.41	0.36	0.43	42.040360	-77.237260	507.6140	Addison
1	ANDE	2025-03-27 15:30:00	-0.4	-1.4	53.1	0.0	0.89	0.0	2.6	5.6	...	2.6	3.9	4.6	0.23	0.13	0.09	42.182270	-74.801390	518.2820	Andes
2	BATA	2025-03-27 15:30:00	3.4	2.6	51.8	0.0	0.00	0.0	5.6	9.5	...	1.6	3.1	4.2	0.28	0.27	0.27	43.019940	-78.135660	276.1200	Batavia
3	BEAC	2025-03-27 15:30:00	5.1	3.9	37.6	0.0	0.00	0.0	1.9	6.2	...	6.4	5.9	6.2	0.40	0.32	0.32	41.528750	-73.945270	90.1598	Beacon
4	BELD	2025-03-27 15:30:00	-0.1	-0.6	51.2	0.0	0.00	0.0	4.4	7.3	...	2.8	3.6	4.3	0.38	0.48	0.40	42.223220	-75.668520	470.3700	Belden
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
122	WFMB	2025-03-27 15:30:00	0.2	-0.4	31.4	0.0	0.26	0.0	0.8	4.8	...	1.6	5.2	3.4	0.29	0.26	0.30	44.393236	-73.858829	614.5990	Whiteface Mountain Base
123	WGAT	2025-03-27 15:30:00	1.5	0.8	36.9	0.0	0.00	0.0	1.5	3.0	...	0.5	1.1	1.4	0.16	0.28	0.10	43.532408	-75.158597	442.9660	Woodgate
124	WHIT	2025-03-27 15:30:00	5.4	4.4	31.9	0.0	0.00	0.0	4.1	6.9	...	2.9	3.2	3.5	0.36	0.50	0.46	43.485073	-73.423071	36.5638	Whitehall
125	WOLC	2025-03-27 15:30:00	3.2	2.5	38.3	0.0	0.00	0.0	2.8	5.0	...	3.0	3.1	4.2	0.19	0.04	0.10	43.228680	-76.842610	121.2190	Wolcott
126	YORK	2025-03-27 15:30:00	4.1	3.3	49.5	0.0	0.00	0.0	4.2	5.5	...	2.9	3.9	4.7	0.25	0.27	0.31	42.855040	-77.847760	177.9420	York

127 rows × 34 columns

Specify a couple of sites using Pandas’ query function.

df.query('station == "LAUR"')

	station	time	temp_2m [degC]	temp_9m [degC]	relative_humidity [percent]	precip_incremental [mm]	precip_local [mm]	precip_max_intensity [mm/min]	avg_wind_speed_prop [m/s]	max_wind_speed_prop [m/s]	...	soil_temp_05cm [degC]	soil_temp_25cm [degC]	soil_temp_50cm [degC]	soil_moisture_05cm [m^3/m^3]	soil_moisture_25cm [m^3/m^3]	soil_moisture_50cm [m^3/m^3]	lat	lon	elevation	name
62	LAUR	2025-03-27 15:30:00	0.1	-0.5	52.0	0.0	0.0	0.0	3.6	6.0	...	1.9	3.1	3.5	0.51	0.48	0.44	42.52526	-75.1369	495.971	Laurens

1 rows × 34 columns

df.query('station == "LKPL"')

	station	time	temp_2m [degC]	temp_9m [degC]	relative_humidity [percent]	precip_incremental [mm]	precip_local [mm]	precip_max_intensity [mm/min]	avg_wind_speed_prop [m/s]	max_wind_speed_prop [m/s]	...	soil_temp_05cm [degC]	soil_temp_25cm [degC]	soil_temp_50cm [degC]	soil_moisture_05cm [m^3/m^3]	soil_moisture_25cm [m^3/m^3]	soil_moisture_50cm [m^3/m^3]	lat	lon	elevation	name
63	LKPL	2025-03-27 15:30:00	-1.4	-1.5	37.9	0.0	0.11	0.0	5.4	7.3	...	0.3	0.9	1.5	0.23	0.28	0.1	44.241615	-73.983454	615.1	Lake Placid

1 rows × 34 columns

Examine the column names.

df.columns

Index(['station', 'time', 'temp_2m [degC]', 'temp_9m [degC]',
       'relative_humidity [percent]', 'precip_incremental [mm]',
       'precip_local [mm]', 'precip_max_intensity [mm/min]',
       'avg_wind_speed_prop [m/s]', 'max_wind_speed_prop [m/s]',
       'wind_speed_stddev_prop [m/s]', 'wind_direction_prop [degrees]',
       'wind_direction_stddev_prop [degrees]', 'avg_wind_speed_sonic [m/s]',
       'max_wind_speed_sonic [m/s]', 'wind_speed_stddev_sonic [m/s]',
       'wind_direction_sonic [degrees]',
       'wind_direction_stddev_sonic [degrees]', 'solar_insolation [W/m^2]',
       'station_pressure [mbar]', 'snow_depth [cm]', 'frozen_soil_05cm [bit]',
       'frozen_soil_25cm [bit]', 'frozen_soil_50cm [bit]',
       'soil_temp_05cm [degC]', 'soil_temp_25cm [degC]',
       'soil_temp_50cm [degC]', 'soil_moisture_05cm [m^3/m^3]',
       'soil_moisture_25cm [m^3/m^3]', 'soil_moisture_50cm [m^3/m^3]', 'lat',
       'lon', 'elevation', 'name'],
      dtype='object')

Create objects pointing to some columns of interest. In Pandas, we can refer to columns by using a “.” in addtion to “[]” in most circumstances, though not if a column name starts with a number, nor if there are spaces in the column name.

stid = df.station
lat = df.lat
lon = df.lon
tmp2 = df['temp_2m [degC]'] # Use brackets due to the presence of a space in the column name
tmp9 = df['temp_9m [degC]']
time = df.time

We will plot one of the variables on a map later on in this notebook. In case we want to go back later and pick a different variable to plot (e.g. 9 m temperature), let’s define a generic object name.

param = tmp9 # replace with tmp9, e.g., if you want to plot 9 m temperature

param

0      1.9
1     -1.4
2      2.6
3      3.9
4     -0.6
      ... 
122   -0.4
123    0.8
124    4.4
125    2.5
126    3.3
Name: temp_9m [degC], Length: 127, dtype: float64

Let’s look at the time object.

time

0      2025-03-27 15:30:00
1      2025-03-27 15:30:00
2      2025-03-27 15:30:00
3      2025-03-27 15:30:00
4      2025-03-27 15:30:00
              ...         
122    2025-03-27 15:30:00
123    2025-03-27 15:30:00
124    2025-03-27 15:30:00
125    2025-03-27 15:30:00
126    2025-03-27 15:30:00
Name: time, Length: 127, dtype: object

The times are the same for all stations, so let’s just pull out one of them, and then create some formatted datetime strings out of it.

timeString = time[0]
timeString

'2025-03-27 15:30:00'

The timeString above could make for a perfectly fine part of an informative title for a map, but let’s create a string of the form “Month Day, Year, HourMin UTC” (e.g., Mar 30, 2021, 0120 UTC )

First, make a datetime object from the string, using strptime with a format that matches that of the string.

timeObj = datetime.strptime(timeString,"%Y-%m-%d %H:%M:%S")
timeObj

datetime.datetime(2025, 3, 27, 15, 30)

Next, make a string object from the datetime object, using strftime with a format that matches what we want to use in the map’s title.

titleString = datetime.strftime(timeObj,"%B %d %Y, %H%M UTC")

titleString

'March 27 2025, 1530 UTC'

Create a scatterplot to show the locations of each NYS Mesonet site using Matplotlib’s scatter method. This method accepts an entire array of lon-lat values.

ax.set_title ('New York State Mesonet Site Locations')
ax.scatter(lon,lat,s=9,c='r',edgecolor='black',alpha=0.75,transform=ccrs.PlateCarree())
# Plot the figure, now with the sites plotted
fig

Did you notice the transform argument? Since we are plotting on a Lambert Conformal-projected map, which uses a Cartesian x-y coordinate system where each point is equally separated in meters, we need to convert, or transform, the lat-lon coordinates into their equivalent coordinates in our chosen projection. We use the transform argument, and assign its value to the coordinate system that our lat-lon array is derived from.

Next, plot the site IDs, using Matplotlib’s text method. This method only accepts a single value for its x and y coordinates, so we need to loop over all the values in the arrays.

Python’s enumerate function:

In a Python list, one can use the enumerate function to set a numerical value to be used as a counter variable (see, e.g… https://realpython.com/python-enumerate/), while at the same time assigning the corresponding list element’s value to another variable. We can use the same technique on a Pandas series object:

for count, value in enumerate(stid):
    print (count, value)

0 ADDI
1 ANDE
2 BATA
3 BEAC
4 BELD
5 BELL
6 BELM
7 BERK
8 BING
9 BKLN
10 BRAN
11 BREW
12 BROC
13 BRON
14 BROO
15 BSPA
16 BUFF
17 BURD
18 BURT
19 CAMD
20 CAPE
21 CHAZ
22 CHES
23 CINC
24 CLAR
25 CLIF
26 CLYM
27 COBL
28 COHO
29 COLD
30 COPA
31 COPE
32 CROG
33 CSQR
34 DELE
35 DEPO
36 DOVE
37 DUAN
38 EAUR
39 EDIN
40 EDWA
41 ELDR
42 ELLE
43 ELMI
44 ESSX
45 FAYE
46 FRED
47 GABR
48 GFAL
49 GFLD
50 GROT
51 GROV
52 HAMM
53 HARP
54 HARR
55 HART
56 HERK
57 HFAL
58 ILAK
59 JOHN
60 JORD
61 KIND
62 LAUR
63 LKPL
64 LOUI
65 MALO
66 MANH
67 MEDI
68 MEDU
69 MORR
70 NBRA
71 NEWC
72 NHUD
73 OLDF
74 OLEA
75 ONTA
76 OPPE
77 OSCE
78 OSWE
79 OTIS
80 OWEG
81 PENN
82 PHIL
83 PISE
84 POTS
85 QUEE
86 RAND
87 RAQU
88 REDF
89 REDH
90 ROXB
91 RUSH
92 SARA
93 SBRI
94 SCHA
95 SCHO
96 SCHU
97 SCIP
98 SHER
99 SOME
100 SOUT
101 SPRA
102 SPRI
103 STAT
104 STEP
105 STON
106 SUFF
107 TANN
108 TICO
109 TULL
110 TUPP
111 TYRO
112 VOOR
113 WALL
114 WALT
115 WANT
116 WARS
117 WARW
118 WATE
119 WBOU
120 WELL
121 WEST
122 WFMB
123 WGAT
124 WHIT
125 WOLC
126 YORK

In the above example, count represents the row index number in the Pandas series, while value is its corresponding value. E.g., the first element’s count is 0, while the value is ADDI. The for loop iterates through all (127) elements of the series.

Now, we can iterate over all 127 sites. For each site, we use Matplotlib’s ax.text function to plot the site ID at its corresponding longitude and latitude. E.g., the first site (ADDI, when count = 0)’s longitude is -77.237260 and its latitude is 42.040360.

for count, value in enumerate(stid):
    ax.text(lon[count],lat[count],value,horizontalalignment='right',transform=ccrs.PlateCarree(),fontsize=7)

fig

Now, let’s attempt to plot the site locations again, but this time we’ll omit the transform argument in ax.scatter and ax.text.

fig = plt.figure(figsize=(11,8.5),dpi=125)
ax = plt.subplot(1,1,1,projection=proj)
ax.set_extent ([lonW,lonE,latS,latN])
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res))
ax.set_title ('New York State Mesonet Site Locations')
ax.scatter(lon,lat,s=9,c='r',edgecolor='black',alpha=0.75)
for count, site in enumerate(stid):
    ax.text(lon[count],lat[count],site,horizontalalignment='right',fontsize=7)

What do you think happened here?

Exercise

Create a new figure which plots the current 2m temperature at the NYSM sites. Include the current date/time to the title.

# %load '/spare11/atm350/common/mar25/cartopy_2.py'
proj = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)
res = '50m' # Use the medium-resolution shapefiles

fig = plt.figure(figsize=(15,10),dpi=125)
ax = plt.subplot(1,1,1,projection=proj)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_title(f"NYS Mesonet 2 meter temperature (°C), {titleString}")
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res));

for count, value in enumerate(param):
    ax.text(lon[count],lat[count],value,horizontalalignment='right',transform=ccrs.PlateCarree(),fontsize=7,color='blue')

Now, plot a different variable … 9 meter temperature. We just need to reset param’s value so it points to the 9-m temperature series; we can repeat the exact same cell as we used above, except we do need to modify the figure title so it reflects what variable we are using.

param = tmp9

fig = plt.figure(figsize=(15,10),dpi=125)
ax = plt.subplot(1,1,1,projection=proj)
ax.set_extent ([lonW,lonE,latS,latN])
ax.set_title(f"NYS Mesonet 9 meter temperature (°C), {titleString}")
ax.add_feature(cfeature.COASTLINE.with_scale(res))
ax.add_feature (cfeature.STATES.with_scale(res));

for count, value in enumerate(param):
    ax.text(lon[count],lat[count],value,horizontalalignment='right',transform=ccrs.PlateCarree(),fontsize=7,color='blue')

We could plot additional variables, and in so doing create a standard surface station plot, but we would find it challenging to align all the variables around the points corresponding to the NYS Mesonet sites. Fortunately, Metpy’s StationPlot method takes care of this, and we will be using it in upcoming notebooks!

Sneak preview: Latest NYSM Surface Map