Geographic Coordinates

Most gridders and processing methods in Verde operate under the assumption that the data coordinates are Cartesian. To process data in geographic (longitude and latitude) coordinates, we must first project them. There are different ways of doing this in Python but most of them rely on the PROJ library. We’ll use pyproj to access PROJ directly and handle the projection operations.

import pyproj
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import verde as vd

With pyproj, we can create functions that will project our coordinates to and from different coordinate systems. For our Baja California bathymetry data, we’ll use a Mercator projection.

data = vd.datasets.fetch_baja_bathymetry()
# We're choosing the latitude of true scale as the mean latitude of our dataset.
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())

The Proj object is a callable (meaning that it behaves like a function) that will take longitude and latitude and return easting and northing coordinates.

# pyproj doesn't play well with Pandas so we need to convert to numpy arrays
proj_coords = projection(data.longitude.values, data.latitude.values)
print(proj_coords)

Out:

(array([-11749316.65303045, -11748999.90777981, -11748682.14077028, ...,
       -11749178.71558258, -11749531.2223938 , -11749882.70744612]), array([2905766.05183735, 2905457.83817751, 2905148.48639437, ...,
       2579612.89349502, 2580017.48641078, 2580422.09116872]))

We can plot our projected coordinates using matplotlib.

plt.figure(figsize=(7, 6))
plt.title("Projected coordinates of bathymetry measurements")
# Plot the bathymetry data locations as black dots
plt.plot(proj_coords[0], proj_coords[1], ".k", markersize=0.5)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")
plt.tight_layout()
plt.show()

Cartesian grids

Now we can use verde.BlockReduce and verde.Spline on our projected coordinates. We’ll specify the desired grid spacing as degrees and convert it to Cartesian using the 1 degree approx. 111 km rule-of-thumb.

spacing = 10 / 60
reducer = vd.BlockReduce(np.median, spacing=spacing * 111e3)
filter_coords, filter_bathy = reducer.filter(proj_coords, data.bathymetry_m)
spline = vd.Spline().fit(filter_coords, filter_bathy)

If we now call verde.Spline.grid we’ll get back a grid evenly spaced in projected Cartesian coordinates.

grid = spline.grid(spacing=spacing * 111e3, data_names=["bathymetry"])
print("Cartesian grid:")
print(grid)

Out:

Cartesian grid:
<xarray.Dataset>
Dimensions:     (easting: 54, northing: 61)
Coordinates:
  * easting     (easting) float64 -1.175e+07 -1.173e+07 ... -1.077e+07
  * northing    (northing) float64 2.074e+06 2.093e+06 ... 3.171e+06 3.19e+06
Data variables:
    bathymetry  (northing, easting) float64 -3.635e+03 -3.727e+03 ... 8.87e+03
Attributes:
    metadata:  Generated by Spline(damping=None, engine='auto',\n    force_co...

We’ll mask our grid using verde.distance_mask to get rid of all the spurious solutions far away from the data points.

grid = vd.distance_mask(proj_coords, maxdist=30e3, grid=grid)

plt.figure(figsize=(7, 6))
plt.title("Gridded bathymetry in Cartesian coordinates")
pc = grid.bathymetry.plot.pcolormesh(cmap="viridis", vmax=0, add_colorbar=False)
plt.colorbar(pc).set_label("bathymetry (m)")
plt.plot(filter_coords[0], filter_coords[1], ".k", markersize=0.5)
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)")
plt.gca().set_aspect("equal")
plt.tight_layout()
plt.show()

Geographic grids

The Cartesian grid that we generated won’t be evenly spaced if we convert the coordinates back to geographic latitude and longitude. Verde gridders allow you to generate an evenly spaced grid in geographic coordinates through the projection argument of the grid method.

By providing a projection function (like our pyproj projection object), Verde will generate coordinates for a regular grid and then pass them through the projection function before predicting data values. This way, you can generate a grid in a coordinate system other than the one you used to fit the spline.

# Get the geographic bounding region of the data
region = vd.get_region((data.longitude, data.latitude))
print("Data region in degrees:", region)

# Specify the region and spacing in degrees and a projection function
grid_geo = spline.grid(
    region=region,
    spacing=spacing,
    projection=projection,
    dims=["latitude", "longitude"],
    data_names=["bathymetry"],
)
print("Geographic grid:")
print(grid_geo)

Out:

Data region in degrees: (245.0, 254.705, 20.0, 29.99131)
Geographic grid:
<xarray.Dataset>
Dimensions:     (latitude: 61, longitude: 59)
Coordinates:
  * longitude   (longitude) float64 245.0 245.2 245.3 ... 254.4 254.5 254.7
  * latitude    (latitude) float64 20.0 20.17 20.33 20.5 ... 29.66 29.82 29.99
Data variables:
    bathymetry  (latitude, longitude) float64 -3.621e+03 ... 9.001e+03
Attributes:
    metadata:  Generated by Spline(damping=None, engine='auto',\n    force_co...

Notice that grid has longitude and latitude coordinates and slightly different number of points than the Cartesian grid.

The verde.distance_mask function also supports the projection argument and will project the coordinates before calculating distances.

grid_geo = vd.distance_mask(
    (data.longitude, data.latitude), maxdist=30e3, grid=grid_geo, projection=projection
)

Now we can use the Cartopy library to plot our geographic grid.

plt.figure(figsize=(7, 6))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Geographic grid of bathymetry")
pc = grid_geo.bathymetry.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), vmax=0, zorder=-1, add_colorbar=False
)
plt.colorbar(pc).set_label("meters")
vd.datasets.setup_baja_bathymetry_map(ax, land=None)
plt.tight_layout()
plt.show()

Total running time of the script: ( 0 minutes 8.756 seconds)