Interpolation in geographic coordinates#

Most interpolators 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 under the hood. We’ll use pyproj to access PROJ directly and handle the projection operations. Verde then offers ways of passing the projection information to use the Cartesian interpolators to generate geographic grids and profiles.

import pandas as pd
import pygmt
import pyproj
import ensaio
import verde as vd

Fetch some data#

To demonstrate how to do this, we’ll use Alps GPS velocity data from ensaio once again:

path_to_data = ensaio.fetch_alps_gps(version=1)
data = pd.read_csv(path_to_data)

We’ll get the data region and add a little padding to it so that our interpolation goes beyond the data bounding box.

region = vd.pad_region(
    vd.get_region((data.longitude, data.latitude)),
    1,  # degree
)

Let’s plot the data on a map with coastlines and country borders to see what it looks like:

fig = pygmt.Figure()
fig.coast(
    region=region, projection="M15c", frame="af",
    land="#eeeeee", borders="1/#666666", area_thresh=1e4,
)
pygmt.makecpt(
    cmap="polar+h",
    series=[data.velocity_up_mmyr.min(), data.velocity_up_mmyr.max()],
)
fig.plot(
    x=data.longitude,
    y=data.latitude,
    fill=data.velocity_up_mmyr,
    style="c0.2c",
    cmap=True,
    pen="0.5p,black",
)
fig.colorbar(frame='af+l"vertical velocity [mm/yr]"')
fig.show()
../_images/geographic_3_0.png

This is much better than just plotting the projected data with no context! Now we can see that the Alps region is moving upward (red dots) and the surrounding regions are moving downward (blue dots).

Projecting data#

For this region, a Mercator projection will do fine since we’re not too close to the poles and the latitude range isn’t very large. We’ll create a projection function and use to project the data coordinates.

projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
coordinates = projection(data.longitude, data.latitude)

We’ve done this before in Making your first grid. We’ll still use the projected coordinates to decimate the data and fit the interpolator. What’s different is how we’ll use the interpolator to make a grid in geographic coordinates.

Fit an interpolator to the Cartesian data#

Use a Cartesian Spline to fit the data, like we did previously.

interpolator = vd.Spline().fit(coordinates, data.velocity_up_mmyr)

Now we can use this interpolator for gridding and predicting a profile.

Make a grid in geographic coordinates#

The interpolator is inherently Cartesian. If we wanted to use to generate a grid in geographic coordinates, we would have to:

  1. Generate grid coordinates on a geographic system.

  2. Project the grid coordinates to Cartesian.

  3. Pass the projected coordinates to the predict method of the interpolator.

  4. Generate an xarray.Dataset with the grid values and the geographic coordinates.

To facilitate this, the grid and profile methods of Verde interpolators take a projection argument. If this is passed, Verde will do the steps above and generate a grid/profile in geographic coordinates. In this case, the region and spacing arguments must be given in geographic coordinates.

grid = interpolator.grid(
    spacing=10 / 60,  # 10 arc-minutes in decimal degrees
    region=region,
    projection=projection,
    dims=("latitude", "longitude"),
    data_names="velocity_up",
)
grid
<xarray.Dataset> Size: 86kB
Dimensions:      (latitude: 76, longitude: 139)
Coordinates:
  * longitude    (longitude) float64 1kB -5.497 -5.329 -5.162 ... 17.42 17.58
  * latitude     (latitude) float64 608B 40.93 41.09 41.26 ... 53.05 53.21 53.38
Data variables:
    velocity_up  (latitude, longitude) float64 85kB -2.309 -2.283 ... -0.2871
Attributes:
    metadata:  Generated by Spline()

Hint

Notice that we set the dims and data_names arguments above. Those control the names of the coordinates and the data variables in the final grid. It’s useful to set those to avoid Verde’s default names, which for this case wouldn’t be appropriate.

Notice that the grid has longitude and latitude coordinates and that they are evenly spaced. We can use this grid to plot a map of the vertical velocity with coastlines and country borders added.

fig = pygmt.Figure()
pygmt.makecpt(
    cmap="polar+h",
    series=[data.velocity_up_mmyr.min(), data.velocity_up_mmyr.max()],
)
fig.grdimage(
    grid.velocity_up,
    cmap=True,
    projection="M15c",
    frame="af",
)
fig.colorbar(frame='af+l"upward velocity (mm/yr)"')
fig.coast(
    shorelines="#333333", borders="1/#666666", area_thresh=1e4,
)
fig.plot(
    x=data.longitude,
    y=data.latitude,
    style="c0.1c",
    fill="#888888",
)
fig.show()
../_images/geographic_7_0.png

Make a profile in geographic coordinates#

Profiles in geographic coordinates would require a similar workflow to grids:

  1. Project the geophysics coordinates of the points to Cartesian.

  2. Generate the profile coordinates using the Cartesian points.

  3. Pass the Cartesian profile coordinates to the predict method of the interpolator.

  4. Convert the projected profile coordinates to geographic with an inverse projection.

Once again, we pass the projection argument to the profile method of the interpolator and let it do the work for us.

profile = interpolator.profile(
    point1=(4, 51),  # longitude, latitude
    point2=(11, 42),
    size=200,  # number of points
    dims=("latitude", "longitude"),
    data_names="velocity_up",
    projection=projection,
)
profile
latitude longitude distance velocity_up
0 51.000000 4.000000 0.000000e+00 -0.020266
1 50.958517 4.035176 5.776034e+03 -0.018964
2 50.916997 4.070352 1.155207e+04 -0.017694
3 50.875440 4.105528 1.732810e+04 -0.016447
4 50.833845 4.140704 2.310414e+04 -0.015208
... ... ... ... ...
195 42.195756 10.859296 1.126327e+06 0.627948
196 42.146874 10.894472 1.132103e+06 0.642185
197 42.097954 10.929648 1.137879e+06 0.656719
198 42.048996 10.964824 1.143655e+06 0.671429
199 42.000000 11.000000 1.149431e+06 0.686197

200 rows × 4 columns

The output comes as a pandas.DataFrame with the longitude and latitude coordinates of the points. The distance is calculated from the projected coordinates and is not a great circle distance.

Lets plot the profile coordinates onto our map and the profile itself to see what it looks like:

fig = pygmt.Figure()
# Plot the grid
pygmt.makecpt(
    cmap="polar+h",
    series=[data.velocity_up_mmyr.min(), data.velocity_up_mmyr.max()],
)
fig.grdimage(grid.velocity_up, cmap=True, projection="M15c", frame="af")
fig.colorbar(frame='af+l"upward velocity (mm/yr)"')
fig.coast(shorelines="#333333", borders="1/#666666", area_thresh=1e4)
fig.plot(
    x=profile.longitude,
    y=profile.latitude,
    fill="#888888",
    style="c0.1c",
)
# Plot the profile above it
fig.shift_origin(yshift="h+1.5c")
fig.plot(
    x=profile.distance,
    y=profile.velocity_up,
    pen="1p",
    projection="X15c/5c",
    frame=["WSne", "xaf+lDistance along profile (m)", "yaf+lUpward velocity (mm/yr)"],
    region=vd.get_region((profile.distance, profile.velocity_up)),
)
fig.show()
../_images/geographic_9_0.png