Raster Rotation and Regridding Concerns

[Hojat-Shirmard]

When reconstructing rasters, GPlately appears to follow this process:

Convert the raster into points (based on a specified resolution)

Rotate the points

Regrid the points using the selected interpolation method

Tasks / Questions:

Check for possible interpolation artifacts (e.g., especially near the poles)

Consider whether the current method (e.g., nearneighbor) is ideal for spherical data

Compare with how GMT handles this — does GMT use a more spherical-aware method?

Feature Request:

Add a default resolution option that automatically matches the input raster (from a NetCDF file or NumPy array) when none is specified. This would simplify usage and help avoid mismatched resolution issues.

[Hojat-Shirmard]

Hi @nickywright @michaelchin @jcannon-gplates

When reconstructing rasters, gplately.raster returns the raster at a default resolution unless the resolution is explicitly defined in the x and y directions. The function signature is:

gplately.raster(
data=None,
plate_reconstruction=None,
extent='global',
realign=False,
resample=None,
resize=None,
time=0.0,
origin=None,
x_dimension_name: str = '',
y_dimension_name: str = '',
data_variable_name: str = '',
**kwargs
)

Also, by default, the interpolation method used is 'linear', and currently, users cannot directly change it to 'nearest', as it's hardcoded internally.

However, we can manually control the interpolation process by following these steps:

1-Convert the raster grid to a set of (lat, lon, value) points.

2-Assign plate IDs using static polygons.

3-Rotate the points

4-Regrid the rotated points using a custom interpolation method—for example, 'nearest' using scipy.interpolate.griddata.

Alternatively, if you'd like, I can try modifying the GPlately source code to support 'nearest' interpolation directly.

[nickywright]

Hi @Hojat-Shirmard,

Which function specifically are you talking about?

Also, I'm confused by some of this (just from the reading of it):

extent='global' - does this mean the raster is assumed to be global? What if it is not? It is then forced to be global? How is it read to be global (I am assuming here (and also from vague memory?) -180 to 180 format, not the 0 to 360 if the grid is in that format - it just doesn't work for half the grid).
Assumptions on extent of rasters is extremely a dangerous assumption in my opinion, especially since you can just read them in...

What is the default resolution? The input?

[Hojat-Shirmard]

Hi @nickywright
Thank you for the feedback!

01-What is the default resolution?
It’s the original resolution of the input raster, unless explicitly overridden.
(This depends on how resize or resample are used:
If resize=(nx, ny) is passed → output is explicitly reshaped to nx × ny.
If resample=(dx, dy) is passed → output is resampled to that grid spacing.
If neither is passed → the output resolution is the same as the input raster's.)

02-Does it assume -180 to 180?
The extent of the input raster is respected unless explicitly overridden.
If the raster is from 0–360 longitude, it stays that way unless realign=True is set, in which case it is realigned to −180 to 180 using realign_grid().
See this function from Gplately in grids.py:

                  def realign_grid(array, lons, lats):
                  mask_lons = lons > 180

                  # realign to -180/180
                  if mask_lons.any():
                      dlon = np.diff(lons).mean()
                      array = np.hstack([array[:, mask_lons], array[:, ~mask_lons]])
                      lons = np.hstack([lons[mask_lons] - 360 - dlon, lons[~mask_lons]])
              
                  if lats[0] > lats[-1]:
                      array = np.flipud(array)
                      lats = lats[::-1]
              
                  return array, lons, lats

03-Does this mean the raster is assumed to be global?
You're right to flag this — the extent='global' keyword does not resample or reproject the grid to a true "global" template. It is just passed downstream and is used more like a label or routing hint.
The extent parameter in the class gplately.Raster(...) only accepts the string 'global' and is not currently implemented to crop or filter spatial bounds.

[Hojat-Shirmard]

We can define a plot function and set the geographic bounds (longitude and latitude) of the map using the extents tuple, which is interpreted with the Plate Carree projection (a basic lat/lon layout).

For example, this line does exactly that:

ax1.set_extent(extents, ccrs.PlateCarree())

Where extents is a tuple of four values:
(min_lon, max_lon, min_lat, max_lat)```
– defining the region to be displayed.

But if we wanna fix this in gplately, we could try modifying gplately.raster—instead of using extent='global', we might use:

extent = ax1.set_extent(extents, ccrs.PlateCarree())

If you’re all on board, I’d be happy to give it a shot!