Properly allow for either pixel or gridline registered grids #476

weiji14 · 2020-06-10T05:27:43Z

Description of proposed changes

xarray assumes pixel registration (aka data values represent centre of the pixel), whereas GMT defaults to using gridline registration (aka data values represent top left corner of pixel). PyGMT currently defaults to gridline registration, following upstream GMT, but users have no way to easily set it to pixel registration.

See http://xarray.pydata.org/en/stable/plotting.html#coordinates and https://docs.generic-mapping-tools.org/latest/cookbook/options.html#grid-registration-the-r-option for context. NOAA also has a good page that talks about grid registration at https://ngdc.noaa.gov/mgg/global/gridregistration.html

This is a major bug, and may cause breaking changes for some users. Users will have to be aware of what type of registration their raster grids' are using, and unfortunately, those using xarray may need to set the registration manually to p (for pixel) for correctness, unless we default to using pixel registration in PyGMT. Alternatively, the salem library has an xarray accessor to represent an xarray grid as a center_grid or corner_grid, as used in weiji14/deepbedmap#150.

Will be good to have some discussion on what's the best plan forward.

Fixes #375? #390? Needed for #451.

Reminders

Run make format and make check to make sure the code follows the style guide.
Add tests for new features or tests that would have caught the bug that you're fixing.
Add new public functions/methods/classes to doc/api/index.rst.
Write detailed docstrings for all functions/methods.
If adding new functionality, add an example to docstrings or tutorials.

Default registration is still "GMT_GRID_NODE_REG" as per upstream, but we might want to override that in some cases.

weiji14 · 2020-06-10T12:56:41Z

Ok, it's getting late in my timezone, but I just want to jot down some notes before bedtime. Please correct me if I'm wrong! But I think we should check Line 105 in the code block below.

pygmt/pygmt/clib/conversion.py

Lines 90 to 113 in 0267dd1

    
           region = [] 
        
           inc = [] 
        
           # Reverse the dims because it is rows, columns ordered. In geographic 
        
           # grids, this would be North-South, East-West. GMT's region and inc are 
        
           # East-West, North-South. 
        
           for dim in grid.dims[::-1]: 
        
               coord = grid.coords[dim].values 
        
               coord_incs = coord[1:] - coord[0:-1] 
        
               coord_inc = coord_incs[0] 
        
               if not np.allclose(coord_incs, coord_inc): 
        
                   raise GMTInvalidInput( 
        
                       "Grid appears to have irregular spacing in the '{}' dimension.".format( 
        
                           dim 
        
                       ) 
        
                   ) 
        
               region.extend([coord.min(), coord.max()]) 
        
               inc.append(coord_inc) 
        
           if any([i < 0 for i in inc]):  # Sort grid when there are negative increments 
        
               inc = [abs(i) for i in inc] 
        
               grid = grid.sortby(variables=list(grid.dims), ascending=True) 
        
           matrix = as_c_contiguous(grid.values[::-1]) 
        
           return matrix, region, inc

Specifically, the code region.extend([coord.min(), coord.max()]) can't surely be valid for both gridline and pixel registered cases. If I'm thinking properly, it is incorrect for the pixel registered case (which xarray uses) by a pixel, the region should be something like coord.min() - inc/2, coord.max() + inc/2).

I've commited a small change at bb7afdc, and have been testing out both GMT_GRID_NODE_REG and GMT_GRID_PIXEL_REG, but they both result in the same figure. Here's a minimal test case I've been playing with.

import pygmt
import xarray as xr

fig = pygmt.Figure()
grid = xr.DataArray(data=[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]])
fig.grdimage(grid=grid, region=[-2, 2, -2, 2], projection="X5c", frame=True)
fig.savefig("pixel_or_node_registration.png")
fig.show()

produces:

This is with the current default gridline settings you can try at https://github.com/GenericMappingTools/try-gmt. Notice how each pixel/square is centred on non-decimal coordinates (i.e. (0,0), (0,1), (0,2), (1,1), etc), and that a 3x3 grid is produced. The same figure is produced if I try using pixel registration, and I feel that it's wrong somehow but my brain can't explain it properly. ~~Will continue tomorrow~~. Edit: Got a bit caught up with other commitments, feel free to jump in anyone if you have any ideas on how to fix this bug.

Trying to use pixel registration as the default when converting an xarray.DataArray to a virtual GMT grid, instead of gridline registration.

weiji14 · 2020-06-19T13:13:23Z

Specifically, the code region.extend([coord.min(), coord.max()]) can't surely be valid for both gridline and pixel registered cases. If I'm thinking properly, it is incorrect for the pixel registered case (which xarray uses) by a pixel, the region should be something like coord.min() - inc/2, coord.max() + inc/2).

Found a way to illustrate this problem: when running grdinfo at https://github.com/GenericMappingTools/try-gmt on the actual grd file vs the virtual earth_relief grid:

import pprint
import pygmt

pprint.pprint(pygmt.grdinfo("@earth_relief_60m"), width=1000)

produces:

('/home/jovyan/.gmt/server/earth_relief_60m.grd: Title: Earth Relief at 01 arc degree\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: Command: grdfilter SRTM15+V2.1.nc -Fg111.2 -D1 -I01d -rp -Gearth/earth_relief/earth_relief_01d_p.grd=ns+s0.5+o0 --IO_NC4_DEFLATION_LEVEL=9 --IO_NC4_CHUNK_SIZE=4096 --PROJ_ELLIPSOID=Sphere\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: Remark: Obtained by Gaussian Cartesian filtering (111.2 km fullwidth) from SRTM15+V2.1.nc [Tozer et al., 2019; http://dx.doi.org/10.1029/2019EA000658]\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: Pixel node registration used [Geographic grid]\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: Grid file format: ns = GMT netCDF format (16-bit integer), CF-1.7\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: x_min: -180 x_max: 180 x_inc: 1 name: longitude n_columns: 360\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: y_min: -90 y_max: 90 y_inc: 1 name: latitude n_rows: 180\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: z_min: -8182 z_max: 5651.5 name: elevation (m)\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: scale_factor: 0.5 add_offset: 0 packed z-range: [-16364,11303]\n'
 '/home/jovyan/.gmt/server/earth_relief_60m.grd: format: netCDF-4 chunk_size: 180,180 shuffle: on deflation_level: 9\n')

whereas

grid = pygmt.datasets.earth_relief.load_earth_relief()
pprint.pprint(pygmt.grdinfo(grid))

produces

(': Title: \n'
 ': Command: \n'
 ': Remark: \n'
 ': Gridline node registration used [Cartesian grid]\n'
 ': Unrecognized grid file format! Probably not a GMT grid\n'
 ': x_min: -179.5 x_max: 179.5 x_inc: 1 name:  n_columns: 360\n'
 ': y_min: -89.5 y_max: 89.5 y_inc: 1 name:  n_rows: 180\n'
 ': z_min: -8182 z_max: 5651.5 name: z\n'
 ': scale_factor: 1 add_offset: 0\n')

Notice how the x_min/x_max of the virtual file isn't at -180/180, ditto with y_min/y_max. Plus it thinks the grid is gridline registered, rather than pixel registered. Also, doesn't seem okay treating the grid as Cartesian instead of Geographic (but that's probably a separate issue with the virtualfile mechanism). Not sure if it's a limitation of the GMT_Create_Data C API? I've pushed a commit using GMT_GRID_PIXEL_REG as the default to see how things turn out. Edit: It runs into a Segmentation Fault 😨

seisman · 2020-06-20T04:00:52Z

There may be a potential bug in GMT. Please follow GenericMappingTools/gmt#3502 and GenericMappingTools/gmt#3503 for upstream feedback.

seisman · 2020-06-20T14:31:11Z

It's confirmed that there was a bug when passing a pixel-registered matrix to GMT C API. It was already fixed in GenericMappingTools/gmt#3503 and merged in to master. I can confirm that the current PR works well. Thus what we need to do is:

decide the default registration for array data
try our best to detect and guess the grid registration (especially for netcdf file loaded as xarray)
let users specify the registration they want to use
wait for the GMT 6.1 release.

seisman · 2020-06-22T20:33:11Z

doesn't seem okay treating the grid as Cartesian instead of Geographic (but that's probably a separate issue with the virtualfile mechanism). Not sure if it's a limitation of the GMT_Create_Data C API?

pygmt/pygmt/clib/session.py

Line 1243 in 86c46f0

family, geometry, mode="GMT_CONTAINER_ONLY", ranges=region, inc=inc

Changing GMT_CONTAINER_ONLY to GMT_CONTAINER_ONLY|GMT_GRID_IS_GEO, then GMT will know it's a geographic grid.

Previously it thought every xarray.DataArray grid was Cartesian! So the load_earth_relief() grid passed to GMT was treated as a Cartesian virtualfile, whereas the actual `@earth_relief` file was properly read as Geographic. Co-authored-by: Dongdong Tian <[email protected]>

seisman · 2020-06-22T22:27:06Z

pygmt/clib/session.py

@@ -1242,7 +1242,7 @@ def virtualfile_from_grid(self, grid, registration="GMT_GRID_PIXEL_REG"):
        gmt_grid = self.create_data(
            family,
            geometry,
-            mode="GMT_CONTAINER_ONLY",
+            mode="GMT_CONTAINER_ONLY|GMT_GRID_IS_GEO",


No. I didn't mean this change. Users can pass both Cartesian and geographic grids to GMT. We need to use GMT_CONTAINER_ONLY or GMT_CONTAINER_ONLY|GMT_GRID_IS_GEO for different grid types.

BTW, GMT_CONTAINER_ONLY is equivalent to GMT_CONTAINER_ONLY|GMT_GRID_IS_CARTESIAN.

Ah sorry, I thought it 'knew' automatically what the grid type was. Should actually get some tests in for this.

BTW, GMT_CONTAINER_ONLY is equivalent to GMT_CONTAINER_ONLY|GMT_GRID_IS_CARTESIAN.

Ok, good to know.

Just for reference, this was added in GMT 5.4. See:

https://docs.generic-mapping-tools.org/dev/changes.html?highlight=gmt_grid_is_geo#api-changes

https://docs.generic-mapping-tools.org/dev/api.html?highlight=gmt_grid_is_geo#exporting-a-data-set

GenericMappingTools/gmt@b6d0af6.

Also need to add it as a valid_modifier.

Should be GMT_IS_OUTPUT instead of GMT_OUTPUT. Also update some docstrings that were missed in the #210 refactor PR.

…ODE_REG" This reverts commit 08edf6f.

For xarray grids read from disk, there is an 'encoding' dictionary, with a 'source' key that gives the path to the file. Running `grdinfo` on that file can give us the grid registration. This works on the earth_relief grids. For grids that are not read from disk, we still default to assuming gridline registration ("GMT_GRID_NODE_REG").

weiji14 · 2020-06-23T04:01:56Z

It's confirmed that there was a bug when passing a pixel-registered matrix to GMT C API. It was already fixed in GenericMappingTools/gmt#3503 and merged in to master. I can confirm that the current PR works well. Thus what we need to do is:
1. decide the default registration for array data

2. try our best to detect and guess the grid registration (especially for netcdf file loaded as xarray)

I've added an automatic gridline/pixel registration detector in 6d4ece4. It simply does grdinfo on the xarray.DataArray's NetCDF source (specifically via grid.encoding["source"], and will work on the earth_relief grids. If that doesn't work (e.g. for grids generated purely in Python), it will fallback to gridline registration by default.

Values go from latitude 89S to 89N. Avoids the `grdimage [INFORMATION]: 359 (of 361) inconsistent grid values at North pole.` message.

…opy" This reverts commit 5a685b2.

weiji14

Implementation-wise, this PR should be done. However, there are cases on Linux when a black image is returned (hard to reproduce consistently, except when running make test), possibly because of a flaky issue with our test suite, or one that might be solved with #517? Not too sure about Windows.

We could merge this in first and try to track down the issue properly (see if it's a problem in PyGMT or upstream with the GMT C API), or wait until we are confident that the bug is isolated. I have kept the test_grdimage_over_dateline.png image small (~10KB) in case we need to change it for any reason.

weiji14 · 2020-07-15T23:12:03Z

pygmt/tests/test_grdimage.py

+@pytest.mark.runfirst
+@pytest.mark.mpl_image_compare
+def test_grdimage_over_dateline(xrgrid):


Running test_grdimage_over_dateline first (before the other tests) makes the tests pass on Linux (see commit 3499047 and bf03e5b). This is a hacky workaround for the flakiness I described at https://github.com/GenericMappingTools/pygmt/pull/476/files#r453429872.

weiji14 · 2020-07-15T23:15:53Z

pygmt/clib/conversion.py

        inc = [abs(i) for i in inc]
        grid = grid.sortby(variables=list(grid.dims), ascending=True)

-    matrix = as_c_contiguous(grid.values[::-1])
+    matrix = as_c_contiguous(grid[::-1].values)


Reverted 5a685b2 at 60c791d, will stick to just grid[::-1].values for now (i.e. flipping the grid left right in xarray rather than in numpy).

seisman · 2020-07-16T03:58:32Z

Now all the tests pass. Does it mean this PR works?

weiji14 · 2020-07-16T04:03:49Z

Now all the tests pass. Does it mean this PR works?

Short answer: Yes (and it has always worked for macOS).

Long answer: See #476 (comment). The tests pass when I run test_grdimage_over_dateline before the other tests. On my Linux box, it works fine if I just call grdimage as normal in a script. But the make test fails consistently for some unknown reason.

seisman · 2020-07-16T04:21:56Z

Makefile

@@ -29,7 +29,8 @@ test:
 	@echo ""
 	@cd $(TESTDIR); python -c "import $(PROJECT); $(PROJECT).show_versions()"
 	@echo ""
-	cd $(TESTDIR); pytest $(PYTEST_ARGS) $(PROJECT)


Please add some comments here, explain why we need two pytest.

seisman · 2020-07-16T04:26:24Z

pygmt/tests/test_grdimage.py

@@ -51,3 +73,19 @@ def test_grdimage_fails():
    fig = Figure()
    with pytest.raises(GMTInvalidInput):
        fig.grdimage(np.arange(20).reshape((4, 5)))
+
+
+@pytest.mark.runfirst


Also add a brief comment here, explain why we need "runfirst".

seisman · 2020-07-16T04:27:41Z

pygmt/tests/test_grdimage.py

+    Ensure no gaps are plotted over the 180 degree international dateline.
+    Specifically checking that `xrgrid.gmt.gtype = 1` sets `GMT_GRID_IS_GEO`,
+    and that `xrgrid.gmt.registration = 0` sets `GMT_GRID_NODE_REG`. Note that
+    there would be a gap over the dateline if a pixel registered grid is used.


Add a link to issue #375?

seisman · 2020-07-16T04:28:39Z

pygmt/tests/test_grdimage.py

+    fig = Figure()
+    assert xrgrid.gmt.registration == 0  # gridline registration
+    xrgrid.gmt.gtype = 1  # geographic coordinate system
+    fig.grdimage(grid=xrgrid, region="g", projection="A0/0/1c", V="d")


Suggested change

fig.grdimage(grid=xrgrid, region="g", projection="A0/0/1c", V="d")

fig.grdimage(grid=xrgrid, region="g", projection="A0/0/1c")

I'll switch it to V="i" for now, so we still get some extra information.

Co-Authored-By: Dongdong Tian <[email protected]>

seisman · 2020-07-16T04:59:15Z

There is a pytest warning:

python3.7/site-packages/_pytest/mark/structures.py:327: PytestUnknownMarkWarning: Unknown pytest.mark.runfirst - is this a typo?  You can register custom marks to avoid this warning - for details, see https://docs.pytest.org/en/latest/mark.html
    PytestUnknownMarkWarning,

weiji14 · 2020-07-16T05:16:14Z

Just for completeness, although this PR does fix #375, it doesn't fix @MarkWieczorek's reported bug at #390. I tried running the example but couldn't get the cylindrical 'Q' projection to plot in xarray (bottom row), and the Mollweide 'W' projection doesn't work too (middle row). Using a NetCDF file works fine though (top row). This results are the same when ran on the code in this PR or the current master branch at f401f85, so it might be an upstream GMT 6.1 issue.


This PR	Master branch at `f401f85`

seisman · 2020-07-16T05:20:50Z

I once tried to debug the issue in #390, it seems the boundary conditions are slightly different for an xarray input and a netCDF input. Mostly likely a GMT bug.

weiji14 · 2020-07-16T05:24:32Z

Yeah, the resampling doesn't work properly on virtualfiles perhaps. Will merge and try to track down the problem another time.

Add registration (r) to common options

ebd3e00

weiji14 added the bug Something isn't working label Jun 10, 2020

Allow virtualfile_from_grid to accept a registration argument

bb7afdc

Default registration is still "GMT_GRID_NODE_REG" as per upstream, but we might want to override that in some cases.

vercel bot temporarily deployed to Preview June 10, 2020 12:31 Inactive

weiji14 mentioned this pull request Jun 19, 2020

Baseline images need updates due to the recent updates of earth relief data #451

Closed

Try using GMT_GRID_PIXEL_REG as default instead of GMT_GRID_NODE_REG

08edf6f

Trying to use pixel registration as the default when converting an xarray.DataArray to a virtual GMT grid, instead of gridline registration.

vercel bot temporarily deployed to Preview June 19, 2020 13:08 Inactive

Merge branch 'master' into grid_registration

8ce734d

vercel bot temporarily deployed to Preview June 22, 2020 03:44 Inactive

vercel bot temporarily deployed to Preview June 22, 2020 22:22 Inactive

seisman reviewed Jun 22, 2020

View reviewed changes

vercel bot temporarily deployed to Preview June 22, 2020 23:02 Inactive

Revert back to using Cartesian mode, but be explicit this time

522ceba

Also need to add it as a valid_modifier.

weiji14 force-pushed the grid_registration branch from 0d61311 to 522ceba Compare June 22, 2020 23:10

vercel bot temporarily deployed to Preview June 22, 2020 23:10 Inactive

Set valid GMT data mode as GMT_IS_OUTPUT

d203cd5

Should be GMT_IS_OUTPUT instead of GMT_OUTPUT. Also update some docstrings that were missed in the #210 refactor PR.

vercel bot temporarily deployed to Preview June 22, 2020 23:14 Inactive

This was referenced Jun 22, 2020

Deciding on default grid registration and grid type for PyGMT #487

Closed

grdimage does not correctly plot xarray DataArrays when they don't contain the redundant longitude at 360 E #375

Closed

Revert "Try using GMT_GRID_PIXEL_REG as default instead of GMT_GRID_N…

90b22a4

…ODE_REG" This reverts commit 08edf6f.

vercel bot temporarily deployed to Preview June 23, 2020 03:16 Inactive

vercel bot temporarily deployed to Preview June 23, 2020 03:40 Inactive

weiji14 added 2 commits July 15, 2020 22:45

Merge branch 'master' into grid_registration

40bc7f9

Reduce image size from 1 inch to 1 cm

d0be1c5

vercel bot temporarily deployed to Preview July 15, 2020 10:49 Inactive

Test with a grid that does not reach the North Pole

46c8471

Values go from latitude 89S to 89N. Avoids the `grdimage [INFORMATION]: 359 (of 361) inconsistent grid values at North pole.` message.

vercel bot temporarily deployed to Preview July 15, 2020 11:37 Inactive

Try just running test_grdimage_over_dateline first before other tests

3499047

vercel bot temporarily deployed to Preview July 15, 2020 21:48 Inactive

Don't use pytest --no-summary which is a pytest 6.0.0rc1 feature

bf03e5b

vercel bot temporarily deployed to Preview July 15, 2020 22:07 Inactive

Revert "Deep copy grid using xarray.DataArray.copy instead of numpy.c…

60c791d

…opy" This reverts commit 5a685b2.

vercel bot temporarily deployed to Preview July 15, 2020 23:14 Inactive

weiji14 marked this pull request as ready for review July 15, 2020 23:28

weiji14 commented Jul 15, 2020

View reviewed changes

seisman reviewed Jul 16, 2020

View reviewed changes

Say why we need to run test_grdimage_over_dateline before other tests

c411a95

Co-Authored-By: Dongdong Tian <[email protected]>

vercel bot temporarily deployed to Preview July 16, 2020 04:42 Inactive

Register the custom runfirst pytest marker

217a7d7

vercel bot temporarily deployed to Preview July 16, 2020 05:03 Inactive

seisman approved these changes Jul 16, 2020

View reviewed changes

weiji14 merged commit 9a4c0a0 into master Jul 16, 2020

weiji14 deleted the grid_registration branch July 16, 2020 05:34

weiji14 mentioned this pull request Aug 10, 2020

Ensure plotting xarray grids with different central meridians work on some projections #560

Merged

5 tasks

weiji14 mentioned this pull request Sep 5, 2020

Revert plotting across dateline workarounds #585

Merged

5 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Properly allow for either pixel or gridline registered grids #476

Properly allow for either pixel or gridline registered grids #476

weiji14 commented Jun 10, 2020 •

edited

Loading

weiji14 commented Jun 10, 2020 •

edited

Loading

weiji14 commented Jun 19, 2020 •

edited

Loading

seisman commented Jun 20, 2020

seisman commented Jun 20, 2020

seisman commented Jun 22, 2020

seisman Jun 22, 2020 •

edited

Loading

weiji14 Jun 22, 2020

weiji14 Jun 22, 2020

weiji14 commented Jun 23, 2020

weiji14 left a comment •

edited

Loading

weiji14 Jul 15, 2020

weiji14 Jul 15, 2020

seisman commented Jul 16, 2020

weiji14 commented Jul 16, 2020

seisman Jul 16, 2020

seisman Jul 16, 2020

seisman Jul 16, 2020

seisman Jul 16, 2020

weiji14 Jul 16, 2020

seisman commented Jul 16, 2020

weiji14 commented Jul 16, 2020

seisman commented Jul 16, 2020

weiji14 commented Jul 16, 2020

	fig.grdimage(grid=xrgrid, region="g", projection="A0/0/1c", V="d")
	fig.grdimage(grid=xrgrid, region="g", projection="A0/0/1c")

Properly allow for either pixel or gridline registered grids #476

Properly allow for either pixel or gridline registered grids #476

Conversation

weiji14 commented Jun 10, 2020 • edited Loading

weiji14 commented Jun 10, 2020 • edited Loading

weiji14 commented Jun 19, 2020 • edited Loading

seisman commented Jun 20, 2020

seisman commented Jun 20, 2020

seisman commented Jun 22, 2020

seisman Jun 22, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

weiji14 commented Jun 23, 2020

weiji14 left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman commented Jul 16, 2020

weiji14 commented Jul 16, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seisman commented Jul 16, 2020

weiji14 commented Jul 16, 2020

seisman commented Jul 16, 2020

weiji14 commented Jul 16, 2020

weiji14 commented Jun 10, 2020 •

edited

Loading

weiji14 commented Jun 10, 2020 •

edited

Loading

weiji14 commented Jun 19, 2020 •

edited

Loading

seisman Jun 22, 2020 •

edited

Loading

weiji14 left a comment •

edited

Loading