Corrupt surface outputs with max_radius and/or gridline registration

nickw · September 27, 2021, 6:24pm

In multiple datasets I am having challenges with corrupt PyGMT surface outputs related to the max_radius and registration options. In both cases the output grid has a diagonal line from lower left to upper right, and shifted and skewed grids that are discontinuous between upper left and lower right halves.

I almost always get the problem when using any setting for max_radius. If I leave out max_radius the problem still persists in some cases but can usually be resolved by switching from registration=g (gridline) to registration=p (pixel) or vice-versa.

In all cases I have applied blockmean to the input data and avoided prime grid cell numbers, and the problem occurs for tiny grids (~4000 data points) and large grids (3.3M data points). I have found that every now and then I get lucky and find a subset of the data where all settings work. The output of blockmean is not corrupted: it is definitely surface.

For grids that have problems, all max_radius settings fail (small or large, cell-based or distance based) and only one of gridline vs pixel will work.

4 examples are shown below for the same input data in standard UTM coordinates on NW-oriented lines:
CORRECT RESULT: max_radius not used (M=None) and pixel registration (registration="p")
n_columns: 304 n_rows: 200

CORRUPT RESULT: max_radius not used (M=None) and gridline registration (registration="g")
n_columns: 305 n_rows: 201

CORRUPT RESULT: max_radius = 5 cells (M="5c") and pixel registration (registration="p")
n_columns: 304 n_rows: 200

CORRUPT RESULT: max_radius = 5 cells (M="5c") and gridline registration (registration="g")
n_columns: 305 n_rows: 201

Any ideas how to avoid these problems and enable max_radius to be used? It looks like the coordinates have been transformed somehow, but then incorrectly transformed back?

maxrjones · September 28, 2021, 8:59pm

I haven’t been able to reproduce this problem using the testing datasets. Would you be able to share a small dataset for which this problem occurs?

A few questions that would help us debug:

What PyGMT version are you using (pygmt.show_versions)?
What operating system are you using?

maxrjones · September 28, 2021, 9:05pm

Nevermind my previous post, I was able to reproduce the problem using the script below.

data = pygmt.datasets.load_sample_bathymetry()
region = pygmt.info(data, spacing=1)
spacing = "5m"
data_median = pygmt.blockmedian(data, region=region, spacing=spacing)
grid = pygmt.surface(data=data_median, region=region, spacing=spacing, M="2c", registration="g")

maxrjones · September 28, 2021, 10:09pm

I think this is a bug. I have opened a GitHub issue over at https://github.com/GenericMappingTools/gmt/issues/5817 and will take a look.

nickw · September 28, 2021, 10:50pm

Awesome! Thanks @maxrjones.

I confirm I see the same results with your test code. Furthermore
grid = pygmt.surface(data=data_median, region=region, spacing=spacing, M=None, registration="g")
also fails.

But this works fine:
grid = pygmt.surface(data=data_median, region=region, spacing=spacing, M=None, registration="p")

maxrjones · September 29, 2021, 3:10pm

@pwessel fixed this bug in https://github.com/GenericMappingTools/gmt/pull/5820. The bug fix will be included in GMT 6.3. For faster use, you could build GMT from source and set up PyGMT to use that version. I think there’s also a way to install the development GMT version with PyGMT using conda with the conda-forge channel, though I have never used that option.

Thanks for reporting the issue!

nickw · September 30, 2021, 10:13pm

That’s great news @maxrjones, thank you. I wasn’t able to build from source, but I see there are conda-forge dev releases available every few weeks. I’ll try one of those when it next updates.

Based on my understanding Paul’s comments about what was fixed there may still be an issue using: max_radius = 5 cells (M="5c" ) and pixel registration (registration="p" ) together? I’ll test it when the new build is available.

weiji14 · October 13, 2021, 2:54am

FYI, we’ve just pushed a new dev release, so you can try using conda install -c conda-forge/label/dev gmt to get the latest patched version of GMT. Let us know if it works or if you need any additional clarification

nickw · October 13, 2021, 5:30pm

Thanks @weiji14!

I was able to test the updates using the same example data @maxrjones used. I confirm that part of the problem was fixed: registration='g' is now recognized correctly.

However there remains a problem when using the combination: M="5c", registration="p":

Using the example code I get the following:

data = pygmt.datasets.load_sample_bathymetry()
region = pygmt.info(data, spacing=1)
spacing = "5m"
data_median = pygmt.blockmedian(data, region=region, spacing=spacing)
grid = pygmt.surface(x=data_median.longitude, y=data_median.latitude, z=data_median.bathymetry,
                     region=region, spacing=spacing, M="5c", registration="p")

But using registration='g' I get the expected result:

data = pygmt.datasets.load_sample_bathymetry()
region = pygmt.info(data, spacing=1)
spacing = "5m"
data_median = pygmt.blockmedian(data, region=region, spacing=spacing)
grid = pygmt.surface(x=data_median.longitude, y=data_median.latitude, z=data_median.bathymetry,
                     region=region, spacing=spacing, M="5c", registration="g")

Both registration methods work as expected when M=None or M is not specified, here with M=None, registration='p':

show_versions() gives:

PyGMT information:
  version: v0.4.1
System information:
  python: 3.9.7 | packaged by conda-forge | (default, Sep 29 2021, 20:33:18)  [Clang 11.1.0 ]
  executable: /Users/nickwilliams/opt/anaconda3/envs/pygmt/bin/python
  machine: macOS-11.6-x86_64-i386-64bit
Dependency information:
  numpy: 1.21.2
  pandas: 1.3.3
  xarray: 0.19.0
  netCDF4: 1.5.7
  packaging: 21.0
  ghostscript: 9.54.0
  gmt: 6.3.0_2f64c67_2021.10.12
GMT library information:
  binary dir: /Users/nickwilliams/opt/anaconda3/envs/pygmt/bin
  cores: 16
  grid layout: rows
  library path: /Users/nickwilliams/opt/anaconda3/envs/pygmt/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/nickwilliams/opt/anaconda3/envs/pygmt/lib/gmt/plugins
  share dir: /Users/nickwilliams/opt/anaconda3/envs/pygmt/share/gmt
  version: 6.3.0

weiji14 · October 13, 2021, 8:49pm

Ok, so maybe the bug wasn’t entirely fixed? @pwessel, could you try to see if this is still a problem (using -M5c and -rp). Equivalent GMT code is as follows:

gmt blockmedian @tut_ship.xyz -R245/255/20/30 -I5m | gmt surface -R245/255/20/30 -I5m -M5c -rp -Goutput.nc
gmt grdimage output.nc -JX5c -png output

I get the same output as the first image in the previous post on GMT 6.3.0_2f64c67_2021.10.12.

Technically, we should be running surface using gridline registration on that sample dataset, and there is a note that surface doesn’t work with pixel registration (see surface — GMT 6.5.0 documentation)… But I’m wondering if we could raise a warning here at least so people don’t get caught out.

Maybe better to run gmt grdsample -T after gmt surface? Note that grdsample isn’t in PyGMT v0.4.1, but will be coming out in PyGMT v0.5.0 (see pygmt.grdsample — PyGMT).

nickw · October 13, 2021, 9:10pm

I saw this, but was confused because pixel registration does seem to work great as long as masking is not used. Or perhaps it is that specific case of pixel registration with masking that the surface docs were trying to avoid? I look forward to trying grdsample in the future.

weiji14 · October 15, 2021, 2:18am

To be honest, I have no idea what is going on I would maybe stick with running surface on gridline registered grids since that seems to be what the algorithm is designed for, and convert to a pixel registered grid afterwards.

Anyways, Paul has fixed the issue in surface must pass -rg when calling grdmask internally by PaulWessel · Pull Request #5868 · GenericMappingTools/gmt · GitHub, which should be out for GMT 6.3.0. The next GMT dev release might take a few more weeks though, so you’ll need to build GMT from source if you want to get the bugfix. In the meantime, if you really insist on using pixel registered grids, you could try using PyGMT grdsample by installing the dev version following Installing — PyGMT

pip install --pre --index-url https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple pygmt

After that, you should be able to resample the gridline registered grid to a pixel registered grid using something like pygmt.grdsample(grid=grid, translate=True). See pygmt.grdsample — PyGMT for more details.