grdsample (pixel to gridline) returns nans

[MarkWieczorek]

Description of the problem

When using the grdsample routine to translate from pixel registration to gridline registration, the output contains nans at the poles. This does not occur when using the GMT command line utilities.

Minimal Complete Verifiable Example

import numpy as np
import pygmt

ppd = 1  # pixels per degree
nlat = 180 * ppd + 1
nlon = 360 * ppd + 1
array = np.ones((nlat, nlon))
spacing = 1 / ppd

x = np.arange(0, 360 + spacing, spacing)
y = np.arange(-90, 90 + spacing, spacing)
xx, yy = np.meshgrid(x, y)

# create gridline registered dataset
gridline = pygmt.xyz2grd(x=xx.flatten(), y=yy.flatten(), z=array.flatten(), spacing=spacing, region=[0, 360, -90, 90], registration='g', coltypes="g")

# convert to pixel registered
pixel = pygmt.grdsample(grid=gridline, translate=True)

# convert back to gridline registered
gridline2 = pygmt.grdsample(grid=pixel, translate=True)

print(gridline)
print(pixel)
print(gridline2)

Full error message

<xarray.DataArray 'z' (lat: 181, lon: 361)>
array([[1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       ...,
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.],
       [1., 1., 1., ..., 1., 1., 1.]], dtype=float32)
Coordinates:
  * lon      (lon) float64 0.0 1.0 2.0 3.0 4.0 ... 356.0 357.0 358.0 359.0 360.0
  * lat      (lat) float64 -90.0 -89.0 -88.0 -87.0 -86.0 ... 87.0 88.0 89.0 90.0
Attributes:
    long_name:     z
    actual_range:  [1. 1.]
<xarray.DataArray 'z' (lat: 180, lon: 360)>
array([[1.008623  , 0.99907345, 1.        , ..., 1.        , 1.        ,
        1.0083693 ],
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        0.99901193],
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        1.        ],
       ...,
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        1.        ],
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        1.        ],
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        1.        ]], dtype=float32)
Coordinates:
  * lon      (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
  * lat      (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
Attributes:
    long_name:     z
    actual_range:  [0.99901193 1.008623  ]
<xarray.DataArray 'z' (lat: 181, lon: 361)>
array([[       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [1.0020888 , 1.0023468 , 0.99943876, ..., 0.99975574, 1.0018516 ,
        1.0020888 ],
       [0.99972796, 0.9997256 , 1.0000663 , ..., 1.0000675 , 0.99908644,
        0.99972796],
       ...,
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        1.        ],
       [1.        , 1.        , 1.        , ..., 1.        , 1.        ,
        1.        ],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan]], dtype=float32)
Coordinates:
  * lon      (lon) float64 0.0 1.0 2.0 3.0 4.0 ... 356.0 357.0 358.0 359.0 360.0
  * lat      (lat) float64 -90.0 -89.0 -88.0 -87.0 -86.0 ... 87.0 88.0 89.0 90.0
Attributes:
    long_name:     z
    actual_range:  [0.99908644 1.00234675]

System information

PyGMT information:
  version: v0.11.0
System information:
  python: 3.12.1 | packaged by conda-forge | (main, Dec 23 2023, 08:01:35) [Clang 16.0.6 ]
  executable: /Users/lunokhod/miniconda3/envs/py312/bin/python
  machine: macOS-14.3-arm64-arm-64bit
Dependency information:
  numpy: 1.26.3
  pandas: 2.1.4
  xarray: 2023.12.0
  netCDF4: 1.6.5
  packaging: 23.2
  contextily: None
  geopandas: None
  ipython: None
  rioxarray: None
  ghostscript: 10.02.1
GMT library information:
  binary version: 6.5.0
  cores: 10
  grid layout: rows
  image layout:
  library path: /Users/lunokhod/miniconda3/envs/py312/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/lunokhod/miniconda3/envs/py312/lib/gmt/plugins
  share dir: /Users/lunokhod/miniconda3/envs/py312/share/gmt
  version: 6.5.0

[MarkWieczorek]

I also note that the GMT commands returns "1"s for all data values for all three grids. In the pixel grid above, the first two rows have values that are not unity.

[seisman]

I believe these are caused by grid padding (extra rows/columns at the grid boundary). Currently, there is very little we can do on the PyGMT side. Hopefully, the issue can be addressed after PR #2398.