grdsample output is incorrect

[seisman]

To reproduce the issue, see below:

In [1]: import pygmt

In [2]: from pygmt.helpers.testing import load_static_earth_relief

In [3]: pygmt.grdsample("@static_earth_relief.nc", region=[-53,-47,-20, -15], spacing=[2, 1])
grdsample [WARNING]: Output sampling interval in x exceeds input interval and may lead to aliasing.
Out[3]:
<xarray.DataArray 'z' (lat: 5, lon: 3)> Size: 120B
array([[460.84375, 482.78125, 848.6875 ],
       [680.46875, 519.09375, 760.84375],
       [867.75   , 579.03125, 804.875  ],
       [551.75   , 666.6875 , 934.09375],
       [411.3125 , 518.4375 , 879.875  ]])
Coordinates:
  * lat      (lat) float64 40B -19.5 -18.5 -17.5 -16.5 -15.5
  * lon      (lon) float64 24B -52.0 -50.0 -48.0
Attributes:
    Conventions:   CF-1.7
    title:
    history:       gmt grdsample @static_earth_relief.nc -G@GMTAPI@-S-O-G-G-G...
    description:
    long_name:     z
    actual_range:  [411.3125  934.09375]

In [4]: grid = load_static_earth_relief()

In [5]: pygmt.grdsample(grid, region=[-53,-47,-20, -15], spacing=[2, 1])
grdsample [WARNING]: Output sampling interval in x exceeds input interval and may lead to aliasing.
Out[5]:
<xarray.DataArray 'z' (lat: 5, lon: 3)> Size: 120B
array([[460.84375, 482.78125, 891.09375],
       [680.46875, 519.09375, 764.9375 ],
       [867.75   , 579.03125, 852.53125],
       [551.75   , 666.6875 , 958.21875],
       [411.3125 , 518.4375 , 931.28125]])
Coordinates:
  * lat      (lat) float64 40B -19.5 -18.5 -17.5 -16.5 -15.5
  * lon      (lon) float64 24B -52.0 -50.0 -48.0
Attributes:
    Conventions:   CF-1.7
    title:
    history:       gmt grdsample @GMTAPI@-S-I-G-M-G-N-000000 -G@GMTAPI@-S-O-G...
    description:
    long_name:     z
    actual_range:  [411.3125  958.21875]