load_earth_relief now returns grid-registrated grids when registration is not specified

[seisman]

Description of the problem

In GMT<=6.3, when users access a remote grid but don't specify the grid registration (e.g., @earth_relief_01d), the pixel-registrated grid is returned by default unless only the gridline-registrated grid is available.

The behavior has changed recently (GenericMappingTools/gmt#6710). The new behavior is well documented in https://docs.generic-mapping-tools.org/dev/datasets/remote-data.html#data-registration.

Optionally, you can append _g or _p to specifically get the gridline-registered or pixel-registered version (if they both exist). If reg is not specified then the behavior depends on whether you are making a plot or processing/extracting a subset of the data:

1. For plots we will return the pixel-registered version unless only the gridline-registered file is available.
2. For grid processing modules we will return the gridline-registered version unless only the pixel-registered file is available. We will also issue a warning since for calculations you should ideally know and specify exactly what you want.

If you do specify a specific registration and that version is not available you will get an error message.

In PyGMT, we provide the load_earth_relief function to load the earth relief dataset into a xarray.DataArray object. When registration is not specified, the function returns a pixel-registrated grid in GMT<=6.3, but will return a gridline-registrated grid in GMT 6.4, thus a few tests will fail.

[maxrjones]

Is the same change needed for load_earth_age? Also, the alternative to providing a note about different behavior between 6.3 and 6.4 is to define the default as gridline in the function signature, which might be cleaner.

[seisman]

Is the same change needed for load_earth_age?

Yes.

Also, the alternative to providing a note about different behavior between 6.3 and 6.4 is to define the default as gridline in the function signature, which might be cleaner.

Sounds a good idea, but for earth_relief dataset, the 15s data only provide the pixel-registrated version and the 03s and 01s data only provide the gridline-registrated version. So registration="gridline" will fail for the 15s data unless we do a careful check.

[willschlitzer]

Is the same change needed for load_earth_age?

Yes.

Also, the alternative to providing a note about different behavior between 6.3 and 6.4 is to define the default as gridline in the function signature, which might be cleaner.

Sounds a good idea, but for earth_relief dataset, the 15s data only provide the pixel-registrated version and the 03s and 01s data only provide the gridline-registrated version. So registration="gridline" will fail for the 15s data unless we do a careful check.

@seisman Would it be best to write in an if statement in load_remote_dataset that checks if gridline is an available registration and, if so, sets it? Something like:

if registration is None:
    if "gridline" in dataset.resolutions[resolution].registrations:
        registration = "gridline"
    else:
        registration = "pixel"

[seisman]

Is the same change needed for load_earth_age?

Yes.

Also, the alternative to providing a note about different behavior between 6.3 and 6.4 is to define the default as gridline in the function signature, which might be cleaner.

Sounds a good idea, but for earth_relief dataset, the 15s data only provide the pixel-registrated version and the 03s and 01s data only provide the gridline-registrated version. So registration="gridline" will fail for the 15s data unless we do a careful check.

@seisman Would it be best to write in an if statement in load_remote_dataset that checks if gridline is an available registration and, if so, sets it? Something like:

if registration is None:
    if "gridline" in dataset.resolutions[resolution].registrations:
        registration = "gridline"
    else:
        registration = "pixel"

Sounds great! Or we can change the order of available registrations in the dataset metadata (i.e., changing Resolution(["pixel", "gridline"], False) to Resolution(["gridline", "pixel"], False), and use:

if registration is None: 
    registration = dataset.resolutions[resolution].registrations[0]

[willschlitzer]

@seisman Is their a way to see what a grid's registration is once it is returned from load_earth_relief? I thought I could just write a test like assert grid.registration == "pixel" but that doesn't appear to be the case.

I added grid.attrs["registration"] = registration to load_remote_dataset, but want to make sure that isn't redundant.

[seisman]

@seisman Is their a way to see what a grid's registration is once it is returned from load_earth_relief? I thought I could just write a test like assert grid.registration == "pixel" but that doesn't appear to be the case.

You can use:

assert grid.gmt.registration == 1 

See https://www.pygmt.org/dev/api/generated/pygmt.GMTDataArrayAccessor.html#pygmt.GMTDataArrayAccessor.