Deciding on default grid registration and grid type for PyGMT

[weiji14]

Description of the problem

PyGMT currently assumes xarray.DataArray grids as gridline registered (GMT_GRID_NODE_REG), and of cartesian type (GMT_GRID_IS_CARTESIAN). This hardcoding has very likely resulted in the issue reported by @MarkWieczorek at #375. Note that this is a non-issue for users passing in NetCDF filenames directly (phew).

Currrently we are working on allowing for either pixel/gridline registered and cartesian/geographic grids at #476, by changing the virtualfile_from_grid mechanism. But we'll need to decide on a good default going forward.

Current (hardcoded) default:

• Gridline registered GMT_GRID_NODE_REG
• Cartesian type GMT_GRID_IS_CARTESIAN

Proposed default:

• Pixel registered GMT_GRID_PIXEL_REG
• Geographic type GMT_GRID_IS_GEO?

For the former, my rationale for changing the default to pixel registration is because xarray does so (see http://xarray.pydata.org/en/stable/plotting.html#coordinates), though there are libraries like salem that allow changing an xarray grid coordinates from corner (gridline) to centre (pixel) based (see salem.Grid.corner_grid).

For the latter (cartesian/geographic), users will probably realize quite quickly that their grid isn't Geographic and set the appropriate Cartesian flag. It will probably be less surprising for those users plotting across the dateline too (e.g. #390).

In any case, I say we keep the current default for now (say for a v0.1.2 release), and introduce any breaking changes in v0.2.0 (also after GMT 6.1.0 comes out).

System information

Please paste the output of python -c "import pygmt; pygmt.show_versions()":

PyGMT information:
  version: v0.1.1+18.gd203cd5f.dirty
System information:
  python: 3.7.6 | packaged by conda-forge | (default, Jun  1 2020, 18:57:50)  [GCC 7.5.0]
  executable: /home/user/miniconda3/envs/pygmt/bin/python
  machine: Linux-5.4.0-37-generic-x86_64-with-debian-bullseye-sid
Dependency information:
  numpy: 1.18.5
  pandas: 1.0.5
  xarray: 0.15.1
  netCDF4: 1.5.3
  packaging: 20.4
  ghostscript: 9.22
  gmt: 6.0.0
GMT library information:
  binary dir: /home/user/miniconda3/envs/pygmt/bin
  cores: 6
  grid layout: rows
  library path: /home/user/miniconda3/envs/pygmt/lib/libgmt.so
  padding: 2
  plugin dir: /home/user/miniconda3/envs/pygmt/lib/gmt/plugins
  share dir: /home/user/miniconda3/envs/pygmt/share/gmt
  version: 6.0.0

[seisman]

Ping @GenericMappingTools/core @GenericMappingTools/python for comments.

[MarkWieczorek]

For the former, my rationale for changing the default to pixel registration is because xarray does so (see http://xarray.pydata.org/en/stable/plotting.html#coordinates),

I would vote to keep the default as gridline registered, in order to be consistent with GMT. To the best of my knowledge, xarray is agnostic about whether the data are gridline or pixel registered, with the sole exception of the plotting routine. In my opinion, xarray should allow you to use either convention, and the fact that they have chosen only one makes it their problem to fix.

[EJFielding]

Yes, the default in GMT has always been gridline registered, so it would make sense to have that as the default in PyGMT.

[weiji14]

Cool, good to hear those thoughts, not complaining really (since it's less work to keep the default :laughing:). Just wanted to voice my concern for users who may not be familiar with GMT (e.g. those coming from climate modelling, remote sensing, etc) but use xarray to handle NetCDF files. I've cross-posted to the forum at https://forum.generic-mapping-tools.org/t/deciding-on-default-grid-registration-and-grid-type-for-pygmt where it will be more visible to other users, and hopefully we'll get a more complete picture across the wider community.

Also, we've recently wrapped an auto-detect gridline/pixel-registration capability into the code at #476, so it should be hopefully be less surprising for xarray users going forward (i.e. they wouldn't have to explicitly choose pixel/gridline to get the right plot). Not to say we shouldn't put some thought into the default behaviour.

Are there any thoughts on the Cartesian/Geographic choice?

[PaulWessel]

I want to add an important point: @weiji14 says that for grid-line registration, "(aka data values represent top left corner of pixel)," this is not actually true. The nodes are always at the center of the "pixel area" that the node represents. It is just that most geophysical/real data are collected at points that align with what we call the gridlines, whereas images collected and possibly translated to grids often use the scheme on the right in that figure. So both systems are widely used for different data sets, there are huge penalties for converting between them (loss of Nyquist and more), and I am sure one's personal preference will be heavily biased by what data set one works with.

[EJFielding]

It is true that GIS and image processing packages all use pixel-registration by default and GMT has been the different one. I guess you have to decide who is the main audience.

I think it makes sense to keep the grids as Cartesian by default. Many times people want to plot something like model output that is a function of distance and depth or some other grid that has nothing to do with a geographic grid. It would be surprising if a simple grid plot suddenly was transformed into a geographic grid.

I have been using QGIS a lot recently, in part because it is much easier to interactively make maps, unlike the shell scripts of GMT. The default registration in QGIS is pixel-registration. It used to have a default geographic lat-long projection for raster data, but in QGIS v. 3, they made it default to a Cartesian grid of unknown project. You have to set the projection to geographic or whatever it is.

Those are my opinions.

[joa-quim]

I was a bit surprise to find that Landsat8 grids (UIn16) come in grid-registration. I guess QGis defaults to pixel because it relies heavily on GDAL for whom everything is pixel-registered only.

[seisman]

@weiji Do you think we should keep gridline/Cartesian as the default and close the issue?

[seisman]

Pinged the wrong guy in my last comment.

Ping @weiji14 again. I think we should keep gridline/Cartesian as the default, which is consistent with GMT's defaults. Closing the issue?

[weiji14]

Yeah, let's keep using GMT's default. Users can change between gridline/pixel or Cartesian/geographic using https://www.pygmt.org/v0.11.0/api/generated/pygmt.GMTDataArrayAccessor.html#pygmt.GMTDataArrayAccessor

P.S. We should look into using gmtlib_get_grdtype as mentioned at https://github.com/GenericMappingTools/pygmt/issues/2194#issuecomment-1368189366

[MarkWieczorek]

Would it be possible to choose geographic automatically if the region spans 0/360/-90/90?

[seisman]

Would it be possible to choose geographic automatically if the region spans 0/360/-90/90?

Yes, it's one of the cases in which we can assume a grid is geographic. GMT has more complicated logic about whether a grid is geographic (see https://github.com/GenericMappingTools/gmt/blob/bc5dacfd24f04687a291727c82a1072e92e91d62/src/gmt_nc.c#L656 if you're interested). We will likely implement the same logic in the next one or two releases (waiting for #2398 first). Please subscribe to #2194 if you want to follow the progress.