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 cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pyproj

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.65303046, -11748999.90777981, -11748682.14077029, ...,
       -11749178.7155826 , -11749531.22239381, -11749882.70744613]), array([2905766.05183735, 2905457.83817751, 2905148.48639437, ...,
       2579612.89349502, 2580017.48641078, 2580422.09116873]))

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

Projected coordinates of bathymetry measurements

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:

/usr/share/miniconda3/envs/test/lib/python3.9/site-packages/verde/base/base_classes.py:463: FutureWarning: The 'spacing', 'shape' and 'region' arguments will be removed in Verde v2.0.0. Please use the 'verde.grid_coordinates' function to define grid coordinates and pass them as the 'coordinates' argument.
  warnings.warn(
Cartesian grid:
<xarray.Dataset>
Dimensions:     (northing: 61, easting: 54)
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()

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)
/usr/share/miniconda3/envs/test/lib/python3.9/site-packages/verde/base/base_classes.py:463: FutureWarning: The 'spacing', 'shape' and 'region' arguments will be removed in Verde v2.0.0. Please use the 'verde.grid_coordinates' function to define grid coordinates and pass them as the 'coordinates' argument.
  warnings.warn(
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()

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.show()

Out:

/usr/share/miniconda3/envs/test/lib/python3.9/site-packages/verde/datasets/sample_data.py:152: UserWarning: All kwargs are being ignored. They are accepted to guarantee backward compatibility.
  warnings.warn(

Profiles¶

For profiles, things are a bit different. The projection is applied to the input points before coordinates are generated. So the profile will be evenly spaced in projected coordinates, not geographic coordinates. This is to avoid issues with calculating distances on a sphere.

The coordinates returned by the profile method will be in geographic coordinates, so projections given to profile must take an inverse argument so we can undo the projection.

The main utility of using a projection with profile is being able to pass in points in geographic coordinates and get coordinates back in that same system (making it easier to plot on a map).

To generate a profile cutting across our bathymetry data, we can use longitude and latitude points taken from the map above).

start = (-114.5, 24.7)
end = (-110, 20.5)
profile = spline.profile(
    point1=start,
    point2=end,
    size=200,
    projection=projection,
    dims=("latitude", "longitude"),
    data_names=["bathymetry"],
)
print(profile)

Out:

      latitude   longitude       distance   bathymetry
  24.700000 -114.500000       0.000000 -4115.540285
  24.679226 -114.477387    3276.548360 -4397.069052
  24.658449 -114.454774    6553.096720 -4766.953633
  24.637668 -114.432161    9829.645080 -5198.947539
  24.616884 -114.409548   13106.193440 -5639.303392
..         ...         ...            ...          ...
20.585669 -110.090452  638926.930221 -3016.467547
20.564257 -110.067839  642203.478581 -2969.351779
20.542841 -110.045226  645480.026941 -2914.069305
20.521422 -110.022613  648756.575301 -2859.557817
20.500000 -110.000000  652033.123661 -2818.101948

[200 rows x 4 columns]

Plot the profile location on our geographic grid from above.

plt.figure(figsize=(7, 6))
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Profile location")
pc = grid_geo.bathymetry.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), vmax=0, zorder=-1, add_colorbar=False
)
plt.colorbar(pc).set_label("meters")
ax.plot(profile.longitude, profile.latitude, "-k", transform=ccrs.PlateCarree())
ax.text(start[0], start[1], "A", transform=ccrs.PlateCarree())
ax.text(end[0], end[1], "B", transform=ccrs.PlateCarree())
vd.datasets.setup_baja_bathymetry_map(ax, land=None)
plt.show()

Out:

/usr/share/miniconda3/envs/test/lib/python3.9/site-packages/verde/datasets/sample_data.py:152: UserWarning: All kwargs are being ignored. They are accepted to guarantee backward compatibility.
  warnings.warn(

And finally plot the profile.

plt.figure(figsize=(8, 3))
ax = plt.axes()
ax.set_title("Profile of bathymetry (A-B)")
ax.plot(profile.distance, profile.bathymetry, "-k")
ax.set_xlabel("Distance (m)")
ax.set_ylabel("Bathymetry (m)")
ax.set_xlim(profile.distance.min(), profile.distance.max())
ax.grid()
plt.tight_layout()
plt.show()

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

Gallery generated by Sphinx-Gallery

Verde v1.7.0

Geographic Coordinates

Contents

Geographic Coordinates¶

Cartesian grids¶

Geographic grids¶

Profiles¶