Introduction to Verde

Unlike xESMF, verde is mainly designed around unstructured grids. Verde has a lot of great examples, we will use an atmospheric one below, but highly reccomend reading though the documents to get a better handle on the API.

https://www.fatiando.org/verde/latest/

Prerequisites

Knowing your way around pandas, xarray, numpy and matplotlib is beneficial. This is not deisgned to be an introduction to any of those packages. We will be using Cartopy for some of the plots as well.

Imports

import pandas as pd
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cf

from appdirs import *

import verde as vd

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

%load_ext watermark
%watermark --iversions

cartopy   : 0.23.0
verde     : 1.8.0
sys       : 3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
matplotlib: 3.8.4
numpy     : 1.26.4
pandas    : 2.2.2

colnames=['remove', 'lon', 'lat', 'date', 'sensor', 'PM25value', 'sensor_type'] 

df = pd.read_csv('../data/airnow_data.csv', names=colnames, header=None)
df = df.drop(columns='remove')
df.head(3)

	lon	lat	date	sensor	PM25value	sensor_type
0	41.783083	-110.537888	2022-09-05T16:00	PM2.5	45	1
1	40.294300	-110.009000	2022-09-05T16:00	PM2.5	40	1
2	40.464722	-109.560830	2022-09-05T16:00	PM2.5	30	1

df.date = pd.to_datetime(df.date)
df.dtypes

lon                   float64
lat                   float64
date           datetime64[ns]
sensor                 object
PM25value               int64
sensor_type             int64
dtype: object

Let’s pick a time that has decent coverage:

df.date.value_counts().head()

date
2022-09-05 16:00:00    36
2022-09-06 07:00:00    36
2022-09-06 20:00:00    36
2022-09-06 19:00:00    36
2022-09-06 18:00:00    36
Name: count, dtype: int64

df2 = df[df.date == '2022-09-06T07:00']
df2.head(3)

	lon	lat	date	sensor	PM25value	sensor_type
540	41.783083	-110.537888	2022-09-06 07:00:00	PM2.5	16	1
541	40.294300	-110.009000	2022-09-06 07:00:00	PM2.5	29	1
542	40.464722	-109.560830	2022-09-06 07:00:00	PM2.5	23	1

Plotting the data

plt.figure(figsize=(8, 5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cf.STATES)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.35, linestyle='--')

gl.xlabels_top = False
gl.ylabels_right = False

gl.xlocator = mticker.FixedLocator([-110, -108,  -106, -104, -102])

plt.scatter(df2.lat, 
            df2.lon, 
            c = df2.PM25value,
            s=50,
            cmap="YlOrBr")

plt.colorbar().set_label("PM 2.5")
plt.show()

/home/runner/miniconda3/envs/gridding-cookbook-dev/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_1_states_provinces_lakes.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Verde Workflows

Some of these workflows might not work super well from a scientific standpoint, this is just to show the mechanics of the verde package.

coordinates = (df2.lon, df2.lat)

Trend Estimation

Same code, different example can be found here: https://www.fatiando.org/verde/latest/tutorials/trends.html

trend = vd.Trend(degree=1).fit(coordinates, df2.PM25value)
print(trend.coef_)

[38.39163129  0.54379964  0.32281663]

/home/runner/miniconda3/envs/gridding-cookbook-dev/lib/python3.10/site-packages/verde/base/least_squares.py:68: FutureWarning: Series.ravel is deprecated. The underlying array is already 1D, so ravel is not necessary.  Use `to_numpy()` for conversion to a numpy array instead.
  regr.fit(jacobian, data.ravel(), sample_weight=weights)

trend_values = trend.predict(coordinates)
residuals = df2.PM25value - trend_values

plt.figure(figsize=(8, 5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cf.STATES)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.35, linestyle='--')

gl.top_labels = False
gl.ylabels_right = False

gl.xlocator = mticker.FixedLocator([-110, -108,  -106, -104, -102])

maxabs = vd.maxabs(residuals)


plt.scatter(df2.lat, 
            df2.lon, 
            c = residuals,
            s=50,
            vmin=-maxabs,
            vmax=maxabs,
            cmap="bwr")

plt.colorbar().set_label("Residuals")
plt.show()

Unlike the Texas example in Verde, hard to see a distinct regional trend.

Data Decimation

spacing = 30 / 60 # Play around with this!

reducer = vd.BlockReduce(reduction=np.median, spacing= spacing)
print(reducer)

BlockReduce(reduction=<function median at 0x7fb078359770>, spacing=0.5)

filter_coords, filter_data = reducer.filter(
    coordinates=(df2.lat, df2.lon), data=df2.PM25value)

/home/runner/miniconda3/envs/gridding-cookbook-dev/lib/python3.10/site-packages/verde/blockreduce.py:177: FutureWarning: Series.ravel is deprecated. The underlying array is already 1D, so ravel is not necessary.  Use `to_numpy()` for conversion to a numpy array instead.
  columns = {"data{}".format(i): comp.ravel() for i, comp in enumerate(data)}
/home/runner/miniconda3/envs/gridding-cookbook-dev/lib/python3.10/site-packages/verde/blockreduce.py:179: FutureWarning: The provided callable <function median at 0x7fb078552ef0> is currently using DataFrameGroupBy.median. In a future version of pandas, the provided callable will be used directly. To keep current behavior pass the string "median" instead.
  blocked = pd.DataFrame(columns).groupby("block").aggregate(reduction)
/home/runner/miniconda3/envs/gridding-cookbook-dev/lib/python3.10/site-packages/verde/blockreduce.py:231: FutureWarning: Series.ravel is deprecated. The underlying array is already 1D, so ravel is not necessary.  Use `to_numpy()` for conversion to a numpy array instead.
  "coordinate{}".format(i): coord.ravel()
/home/runner/miniconda3/envs/gridding-cookbook-dev/lib/python3.10/site-packages/verde/blockreduce.py:236: FutureWarning: The provided callable <function median at 0x7fb078552ef0> is currently using DataFrameGroupBy.median. In a future version of pandas, the provided callable will be used directly. To keep current behavior pass the string "median" instead.
  grouped = table.groupby("block").aggregate(self.reduction)

Sanity check

np.shape(filter_coords)[1] == np.size(filter_data)

True

Plotting the decimated datapoints

plt.figure(figsize=(6, 4))
plt.plot(*filter_coords, ".k",  markersize=15, label='Data points')
plt.legend()
plt.show()

Verde Spline

This turtorial runs though this entire workflow, but with a bathymetery dataset: https://www.fatiando.org/verde/latest/tutorials/projections.html

spline = vd.Spline().fit(filter_coords, filter_data)

region = vd.get_region(filter_coords)
grid_coords = vd.grid_coordinates(region, spacing=spacing)
grid = spline.grid(coordinates=grid_coords, data_names="PM25")
print(grid)

<xarray.Dataset> Size: 2kB
Dimensions:   (northing: 14, easting: 18)
Coordinates:
  * easting   (easting) float64 144B -110.5 -110.0 -109.5 ... -102.4 -101.9
  * northing  (northing) float64 112B 35.14 35.66 36.18 ... 40.84 41.36 41.88
Data variables:
    PM25      (northing, easting) float64 2kB 74.69 69.09 63.62 ... 11.73 6.292
Attributes:
    metadata:  Generated by Spline(mindist=0)

distance_mask = 1 #this might not make physical sense
grid = vd.distance_mask(filter_coords, maxdist=distance_mask, grid=grid)

Let’s make a plot with the gridded and masked output, also shows where the datapoints are:

plt.figure(figsize=(8, 5))

pc = grid.PM25.plot.pcolormesh(cmap="YlOrBr", add_colorbar=False)
plt.colorbar(pc)
plt.plot(filter_coords[0], filter_coords[1], ".k", markersize=15, label='Data Points')

plt.title('Gridded PM2.5 Data, Distance Filter: '+str(distance_mask))

plt.xlabel("West", size=15)
plt.ylabel("North", size=15)
plt.legend()
plt.gca().set_aspect("equal")
plt.tight_layout()
plt.show()

Xarray

Making the grid into an xarray dataset is 1 line of code!

ds = grid.PM25.to_dataset()
ds

<xarray.Dataset> Size: 2kB
Dimensions:   (easting: 18, northing: 14)
Coordinates:
  * easting   (easting) float64 144B -110.5 -110.0 -109.5 ... -102.4 -101.9
  * northing  (northing) float64 112B 35.14 35.66 36.18 ... 40.84 41.36 41.88
Data variables:
    PM25      (northing, easting) float64 2kB nan nan nan nan ... 17.54 nan nan

xarray.Dataset

• Dimensions:
  • easting: 18
  • northing: 14
• Coordinates: (2)
  • easting
    (easting)
    float64
    -110.5 -110.0 ... -102.4 -101.9

    array([-110.537888, -110.030265, -109.522642, -109.01502 , -108.507397,
           -107.999774, -107.492151, -106.984528, -106.476905, -105.969283,
           -105.46166 , -104.954037, -104.446414, -103.938791, -103.431168,
           -102.923546, -102.415923, -101.9083  ])

  • northing
    (northing)
    float64
    35.14 35.66 36.18 ... 41.36 41.88

    array([35.1441  , 35.661904, 36.179709, 36.697513, 37.215317, 37.733122,
           38.250926, 38.76873 , 39.286534, 39.804339, 40.322143, 40.839947,
           41.357752, 41.875556])

• Data variables: (1)
  • PM25
    (northing, easting)
    float64
    nan nan nan nan ... 17.54 nan nan

    metadata :

    Generated by Spline(mindist=0)

    array([[        nan,         nan,         nan,         nan,         nan,
                    nan, 44.59131046, 39.577399  , 30.9290346 , 19.37611031,
            15.19816161,         nan,         nan,         nan,         nan,
                    nan, 25.41683029, 26.43790005],
           [        nan,         nan,         nan,         nan,         nan,
                    nan,         nan, 26.17680332, 15.9536217 ,  6.36852569,
             6.78126555,  9.51742429,         nan,         nan,         nan,
                    nan, 21.73404816, 22.63338388],
           [        nan,         nan,         nan,         nan,         nan,
            28.83239662, 23.08102862,         nan,  9.79006308,  5.30588243,
             6.1206478 ,  8.59800223,         nan,         nan,         nan,
                    nan,         nan, 19.53628538],
           [        nan,         nan,         nan,         nan, 28.57511624,
            22.36068765, 17.44647738, 13.61649596, 10.31282478,  8.86834165,
             9.71426225, 11.30711447,         nan,         nan,         nan,
                    nan,         nan,         nan],
           [        nan,         nan,         nan,         nan, 25.77930382,
            20.34310532, 16.98000413, 15.49218454,         nan, 14.20224356,
            14.85170643,         nan,         nan,         nan,         nan,
                    nan,         nan,         nan],
    ...
           [36.66061715, 34.20234403, 33.18881741, 36.55249975, 40.26405683,
            42.28480502, 44.22890483, 44.98405891, 42.66580598, 37.17976218,
            29.07744676, 23.71653844, 18.99329113, 15.95589749,         nan,
                    nan,         nan,         nan],
           [31.82417211, 28.7639875 , 25.62738039, 31.31531917, 37.16895659,
            40.01590078,         nan,         nan,         nan, 41.2752953 ,
            39.3087373 , 29.65399235, 17.80505315, 16.26805185,         nan,
                    nan,         nan,         nan],
           [25.60939056, 21.23512155, 18.60734344, 22.6149659 , 28.79845456,
                    nan,         nan,         nan,         nan, 41.18429168,
            43.17136022, 40.17053881, 29.23161609, 22.74468372,         nan,
                    nan,         nan,         nan],
           [19.99437318, 15.64671639, 13.05991168, 15.09135734, 20.54979891,
                    nan,         nan,         nan, 35.53659168, 37.03471545,
            39.94845054, 43.97239772, 38.47274977, 29.53796423, 22.51611929,
            16.90156799,         nan,         nan],
           [15.08874012, 11.13853538,  8.65034358,  9.72304681, 14.06563031,
                    nan,         nan,         nan,         nan, 33.19423011,
            36.73937806, 39.70671614, 37.41064739, 30.95188451, 23.921356  ,
            17.53679109,         nan,         nan]])

• Indexes: (2)
  • easting
    PandasIndex

    PandasIndex(Index([        -110.537888, -110.03026517647058, -109.52264235294118,
           -109.01501952941176, -108.50739670588234, -107.99977388235294,
           -107.49215105882352,  -106.9845282352941,  -106.4769054117647,
           -105.96928258823529, -105.46165976470587, -104.95403694117647,
           -104.44641411764705, -103.93879129411765, -103.43116847058823,
           -102.92354564705882, -102.41592282352941,           -101.9083],
          dtype='float64', name='easting'))

  • northing
    PandasIndex

    PandasIndex(Index([           35.1441,  35.66190430769231,  36.17970861538462,
            36.69751292307693,  37.21531723076923,  37.73312153846154,
            38.25092584615385,  38.76873015384616, 39.286534461538466,
           39.804338769230775,  40.32214307692308, 40.839947384615385,
           41.357751692307694,          41.875556],
          dtype='float64', name='northing'))

• Attributes: (0)

Exporting this as a netCDF:

ds.to_netcdf('test.nc')

Summary

Verde and xESMF are two different gridding packages, with two different aims. They both work well with xarray!