Spatial Temperature Modeling with SPDE

2. Coordinate Reference System

For spatial modeling, we need projected coordinates (kilometers) rather than geographic coordinates (degrees):

# Coordinate Reference Systems
CRS_LONGLAT = "EPSG:4326"  # Standard WGS84 (degrees)
CRS_KM = "+proj=utm +zone=40 +ellps=intl +towgs84=-143,-236,7,0,0,0,0 +units=km +no_defs +type=crs"

# Study region
COUNTRIES = ["Lebanon", "Syria", "Jordan", "Saudi Arabia"]

This is essential because:

• The SPDE model assumes distances are in consistent units
• One degree of longitude varies in distance depending on latitude
• UTM Zone 40 is appropriate for the Middle East region

3. Load Temperature Data

We load the monthly temperature data from Part 1:

# Load temperature data from Part 1
tavg = pd.read_csv('tavg_month.csv')
print(f"Temperature data loaded: {tavg.shape[0]} stations")

# First 2 rows:
#        station  latitude  longitude  elevation  X202401  X202402  ...
# 0  BA000041150   26.2670      50.65        2.0    20.71    19.93  ...
# 1  IS000004640   32.9667      35.50      934.0     9.05     8.74  ...

Temperature data loaded: 55 stations

Load country boundaries and project to the km-based CRS:

# Load country boundaries
NE_URL = "https://naciscdn.org/naturalearth/50m/cultural/ne_50m_admin_0_countries.zip"
worldMapll = gpd.read_file(NE_URL)
map_ll = worldMapll[worldMapll["ADMIN"].isin(COUNTRIES)].copy()
map2km = map_ll.to_crs(CRS_KM)

# Create station locations in km CRS
geometry = [Point(lon, lat) for lon, lat in zip(tavg["longitude"], tavg["latitude"])]
tavg_gdf = gpd.GeoDataFrame(tavg, geometry=geometry, crs=CRS_LONGLAT)
tavg_km = tavg_gdf.to_crs(CRS_KM)
coords_km = np.column_stack([tavg_km.geometry.x, tavg_km.geometry.y])

# Plot stations: color = temperature, size = elevation
fig, ax = plt.subplots(figsize=(10, 8))
map2km.boundary.plot(ax=ax, color='black', linewidth=1.5)
sc = ax.scatter(
    coords_km[:, 0], coords_km[:, 1],
    c=tavg['X202401'], cmap='RdYlBu_r',
    s=tavg['elevation'] / 20 + 20,
    edgecolors='black', linewidth=0.5, zorder=5
)
plt.colorbar(sc, ax=ax, shrink=0.7, label='January 2024 mean temperature (°C)')
ax.set_title('Weather stations: color = temperature, size = elevation')
ax.set_xlabel('Easting (km)')
ax.set_ylabel('Northing (km)')
ax.set_aspect('equal')
plt.tight_layout()
plt.show()

5. Create the Triangular Mesh

The SPDE approach requires a triangular mesh that covers the study region. The mesh serves as the computational domain for the spatial random field.

Mesh Construction Steps

1. Buffer the boundary by 100 km (to avoid edge effects)
2. Generate hexagonal lattice points (edge length = 50 km)
3. Build triangular mesh with specified resolution

# Create unified boundary and buffer by 100 km
boundary = unary_union(map2km.geometry)
boundary_buffered = boundary.buffer(100)

# Extract boundary coordinates for hexagonal lattice
if boundary_buffered.geom_type == 'MultiPolygon':
    largest = max(boundary_buffered.geoms, key=lambda p: p.area)
    boundary_coords = np.array(largest.exterior.coords)
else:
    boundary_coords = np.array(boundary_buffered.exterior.coords)

# Create hexagonal lattice points
hexpoints = fm_hexagon_lattice(boundary_coords, edge_len=50)
print(f"Hexagonal lattice: {len(hexpoints)} points")

# Build the triangular mesh
mesh = fm_mesh_2d(
    loc=hexpoints,
    offset=500,     # Extend 500 km beyond boundary
    max_edge=200    # Maximum triangle edge: 200 km
)

print(f"Mesh: {mesh.n} vertices, {mesh.n_triangle} triangles")

Mesh: 1,792 vertices, 3,499 triangles

# Plot mesh using built-in method
fig, ax = mesh.plot(locs=coords_km, title='SPDE Mesh', show=False)
map2km.boundary.plot(ax=ax, color='black', linewidth=2)
plt.tight_layout()
plt.show()

9. Model Fit

# Plot fitted vs observed
fitted_mean = sfit.summary_fitted_values['mean'].values
observed = tavg['X202401'].values
valid = ~np.isnan(observed)
r_sq = np.corrcoef(fitted_mean[valid], observed[valid])[0, 1] ** 2
print(f"R-squared: {r_sq:.4f}")

R-squared: 0.9788

The $R^2 = 0.9788$ indicates the model (intercept + elevation + spatial random effect) explains nearly all of the observed temperature variation.

fig, ax = plt.subplots(figsize=(7, 7))
ax.scatter(fitted_mean[valid], observed[valid], marker='o', s=60,
           c='#3498db', alpha=0.7)
lims = [min(fitted_mean[valid].min(), observed[valid].min()) - 1,
        max(fitted_mean[valid].max(), observed[valid].max()) + 1]
ax.plot(lims, lims, 'r--', linewidth=2, label='Perfect fit')
ax.set_xlabel('Fitted (posterior mean, °C)')
ax.set_ylabel('Observed (°C)')
ax.set_title(f'Model fit: R² = {r_sq:.4f}')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

10. Create the Grid Projector

To visualize spatial fields on a regular grid, we create a grid projector that maps from mesh vertices to a 200 x 200 pixel grid covering the study region:

# Create grid projector for visualization
bb = map2km.total_bounds  # [minx, miny, maxx, maxy]
bb0_x = (bb[0], bb[2])
bb0_y = (bb[1], bb[3])

projGrid = spde_grid_projector(
    mesh=mesh,
    xlim=bb0_x,
    ylim=bb0_y,
    dims=(200, 200)
)
print(f"Grid projector: {projGrid.dims[0]}x{projGrid.dims[1]} = {projGrid.n_grid} points")
print(f"X range: {bb0_x[0]:.0f} to {bb0_x[1]:.0f} km")
print(f"Y range: {bb0_y[0]:.0f} to {bb0_y[1]:.0f} km")

Grid projector: 200x200 = 40000 points
X range: -1730 to 360 km
Y range: 1864 to 4272 km

11. Visualize the Spatial Effect

The spatial random effect $u(s)$ captures geographic patterns not explained by elevation:

# Find grid points outside the country boundary
bnd = unary_union(map2km.geometry)
sfGrid = gpd.GeoDataFrame(
    geometry=gpd.points_from_xy(projGrid.lattice[:, 0], projGrid.lattice[:, 1]),
    crs=CRS_KM
)
igrid_out = ~sfGrid.within(bnd)

# Project spatial effect to grid
spatial_mean = sfit.summary_random['spatial']['mean'].values
gridSmean = projGrid.project(spatial_mean, mask_outside=False)
gridSmean[igrid_out.values.reshape(projGrid.dims[1], projGrid.dims[0])] = np.nan

# Plot spatial effect
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(
    gridSmean, origin='lower',
    extent=[bb0_x[0], bb0_x[1], bb0_y[0], bb0_y[1]],
    cmap='RdBu_r', aspect='equal'
)
map2km.boundary.plot(ax=ax, color='black', linewidth=1.5)
plt.colorbar(im, ax=ax, shrink=0.7, label='Spatial effect (°C)')
ax.set_title('Spatial random effect (posterior mean)')
ax.set_xlabel('Easting (km)')
ax.set_ylabel('Northing (km)')
plt.tight_layout()
plt.show()

What This Tells Us

• Positive values (red): Warmer than expected based on elevation alone
• Negative values (blue): Cooler than expected

This captures regional climate patterns like coastal effects, continental effects, and local topographic effects.

12. Load Elevation Data (ETOPO)

To make predictions across the entire region, we need elevation at every grid point. We use ETOPO 2022 from NOAA via OpenDAP:

import xarray as xr

# Transform grid points back to lon/lat for ETOPO lookup
sfGrid_ll = sfGrid.to_crs(CRS_LONGLAT)
glocs = np.column_stack([sfGrid_ll.geometry.x, sfGrid_ll.geometry.y])

# Get bounding box in lon/lat (with some padding)
lon_min, lon_max = glocs[:, 0].min() - 1, glocs[:, 0].max() + 1
lat_min, lat_max = glocs[:, 1].min() - 1, glocs[:, 1].max() + 1

# Load ETOPO elevation data via OpenDAP
ETOPO_URL = "https://www.ngdc.noaa.gov/thredds/dodsC/global/ETOPO2022/60s/60s_surface_elev_netcdf/ETOPO_2022_v1_60s_N90W180_surface.nc"

ds = xr.open_dataset(ETOPO_URL, engine='netcdf4')
etopo = ds['z'].sel(
    lon=slice(lon_min, lon_max),
    lat=slice(lat_min, lat_max)
)

# Interpolate elevation at each grid point location
etopoll = etopo.interp(
    lon=xr.DataArray(glocs[:, 0], dims='points'),
    lat=xr.DataArray(glocs[:, 1], dims='points'),
    method='nearest'
).values

# Mask grid points outside the boundary
etopoll[igrid_out.values] = np.nan

ETOPO provides:

• Global coverage (land elevation + ocean depth)
• 1 arc-minute resolution (~1.8 km)
• Streaming access, no download required

Summary

References

• Lindgren, F., Rue, H., and Lindstrom, J. (2011). "An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach." Journal of the Royal Statistical Society: Series B, 73(4), 423-498.
• Fuglstad, G.-A., Simpson, D., Lindgren, F., and Rue, H. (2019). "Constructing Priors that Penalize the Complexity of Gaussian Random Fields." Journal of the American Statistical Association, 114(525), 445-452.
• GHCN-Daily: Menne, M.J., et al. (2012). "Global Historical Climatology Network - Daily." NOAA National Climatic Data Center.