deprecated `pygmt.load_dataarray()` inconsistent with `xr.load_datarray`

[mdtanker]

Description of the problem

I'm not sure if this is a bug or a change to PyGMT which I missed. Since the pygmt.load_dataarray function is deprecated, I'm using the suggested alternative, xr.load_datarray(..., engine="gmt", raster_kind='grid'), however, this gives an incorrect grid registration.

As you can see below, a grid with a registration of 0 (gridline), when converted to pixle registered with grdedit and then opened again, gives the correct pixel registration (1) when using pygmt.load_datarray, but an incorrect gridline registration (0) when opened with xr.load_dataarray.

Minimal Complete Verifiable Example

import pygmt
import xarray as xr

grid = pygmt.datasets.load_earth_relief(
    resolution="30s", region=[144.5, 145.5, 13, 14.5], registration="gridline"
)
print(grid.gmt.registration)
---> 0

with pygmt.clib.Session() as ses:
    with ses.virtualfile_from_grid(grid) as f_in:
        with pygmt.helpers.GMTTempFile(suffix=".nc") as tmpfile:
            args = f"{f_in} -T -G{tmpfile.name}"
            ses.call_module("grdedit", args)
            f_out = pygmt.load_dataarray(tmpfile.name)
print(f_out.gmt.registration)
---> 1 # this is correct

with pygmt.clib.Session() as ses:
    with ses.virtualfile_from_grid(grid) as f_in:
        with pygmt.helpers.GMTTempFile(suffix=".nc") as tmpfile:
            args = f"{f_in} -T -G{tmpfile.name}"
            ses.call_module("grdedit", args)
            f_out = xr.load_dataarray(tmpfile.name, engine='gmt', raster_kind='grid')
print(f_out.gmt.registration)
---> 0  # this is incorrect

Full error message

System information

PyGMT information:
  version: v0.16.0
System information:
  python: 3.12.11 | packaged by conda-forge | (main, Jun  4 2025, 14:45:31) [GCC 13.3.0]
  executable: /home/sungw937/miniforge3/envs/polartoolkit/bin/python
  machine: Linux-5.15.167.4-microsoft-standard-WSL2-x86_64-with-glibc2.39
Dependency information:
  numpy: 2.2.6
  pandas: 2.3.1
  xarray: 2025.7.1
  packaging: 25.0
  contextily: None
  geopandas: 1.1.1
  IPython: 9.4.0
  pyarrow: 20.0.0
  rioxarray: 0.19.0
  gdal: None
  ghostscript: 10.04.0
GMT library information:
  version: 6.5.0
  padding: 2
  share dir: /home/sungw937/miniforge3/envs/polartoolkit/share/gmt
  plugin dir: /home/sungw937/miniforge3/envs/polartoolkit/lib/gmt/plugins
  library path: /home/sungw937/miniforge3/envs/polartoolkit/lib/libgmt.so
  cores: 8
  grid layout: rows
  image layout: 
  binary version: 6.5.0

[weiji14]

Hi @mdtanker! Hope you're settling well at your new postdoc position 😄

I can reproduce the behaviour you posted. What is happening is that we are now dynamically inferring the grid registration from the source encoding file, so if you are using a temporary file that is deleted, there is no way for PyGMT to retrieve the grid registation info, and it fallbacks to the default of 0.

To illustrate, if you put the print statement inside the with-block:

with pygmt.clib.Session() as ses:
    with ses.virtualfile_from_grid(grid) as f_in:
        with pygmt.helpers.GMTTempFile(suffix=".nc") as tmpfile:
            args = f"{f_in} -T -G{tmpfile.name}"
            ses.call_module("grdedit", args)
            f_out = xr.load_dataarray(tmpfile.name, engine='gmt', raster_kind='grid')
            print(f_out.gmt.registration)
---> 1  # this is correct

whereas if you put it outside, where the tempfile is deleted already:

with pygmt.clib.Session() as ses:
    with ses.virtualfile_from_grid(grid) as f_in:
        with pygmt.helpers.GMTTempFile(suffix=".nc") as tmpfile:
            args = f"{f_in} -T -G{tmpfile.name}"
            ses.call_module("grdedit", args)
            f_out = xr.load_dataarray(tmpfile.name, engine='gmt', raster_kind='grid')
print(f_out.gmt.registration)
---> 0  # this is incorrect

alternative workaround, if you insist on putting the print() statement outside the with-block, is to cache the gmt accessor info like this:

with pygmt.clib.Session() as ses:
    with ses.virtualfile_from_grid(grid) as f_in:
        with pygmt.helpers.GMTTempFile(suffix=".nc") as tmpfile:
            args = f"{f_in} -T -G{tmpfile.name}"
            ses.call_module("grdedit", args)
            f_out = xr.load_dataarray(tmpfile.name, engine='gmt', raster_kind='grid')
            _ = f_out.gmt # cache gmt accessor info
print(f_out.gmt.registration)
---> 1  # this is correct

Looking at the code here:

https://github.com/GenericMappingTools/pygmt/blob/eddc912cbfc918a709ce96efb0ac77d2b3e44179/pygmt/xarray/backend.py#L148-L153

It seems like we're not caching the gmt accessor values (registation and gtype) for non-virtual files. I took out the _ = raster.gmt line at https://github.com/GenericMappingTools/pygmt/pull/3932#discussion_r2067584046, because I thought lib.virtualfile_to_raster would cache the values, but that works only for virtual files, not non-virtual files like here. So I think yes, there is a bug here, that only happens when a) you passed in the grid as a NetCDF (non-virtualfile), and b) the file is deleted before the accessor info is retrieved. What do you think @seisman?

[seisman]

I'm still trying to understand what's exactly happening.

@mdtanker For your case, I would recommend avoid using temporary files if possible. The following script works as expected:

import pygmt
import xarray as xr

grid = pygmt.datasets.load_earth_relief(
    resolution="30s", region=[144.5, 145.5, 13, 14.5], registration="gridline"
)

with pygmt.clib.Session() as ses:
    with ses.virtualfile_in(data=grid) as vingrd, ses.virtualfile_out(kind="grid") as voutgrd:
        args = f"{vingrd} -T -G{voutgrd}"
        ses.call_module("grdedit", args)
        f_out =  ses.virtualfile_to_raster(vfname=voutgrd, kind="grid")
print(f_out.gmt.registration)

[mdtanker]

Thank you both for your fast responses! I'm using @seisman's solution and it's working well.

Thanks @weiji14! Yeah, feels good to be able to settle down and not be on the job hunt 😆

[yvonnefroehlich]

Thanks @weiji14! Yeah, feels good to be able to settle down and not be on the job hunt 😆

Great news @mdtanker ! I really hope you will enjoy your time in Kiel 🚀! I worked with Jörg and Wolfgang during my Bachelor thesis 🙂 (and started learning and using GMT).

[mdtanker]

@yvonnefroehlich Thanks, I'm definitely enjoying it so far! Small world 😆