Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Geographic xarray.DataArray grid is plotted as Cartesian grid #1172

Closed
seisman opened this issue Apr 4, 2021 · 3 comments · Fixed by #1224
Closed

Geographic xarray.DataArray grid is plotted as Cartesian grid #1172

seisman opened this issue Apr 4, 2021 · 3 comments · Fixed by #1224
Labels
bug Something isn't working upstream Bug or missing feature of upstream core GMT
Milestone

Comments

@seisman
Copy link
Member

seisman commented Apr 4, 2021

Description of the problem

Originally found in #1154 (comment)

Full code that generated the error

import pygmt

pygmt.show_versions()

xrgrid = pygmt.grdcut("@earth_relief_01d", region=[-116, -109, -47, -44])
print(xrgrid.gmt.gtype)

fig = pygmt.Figure()
fig.grdview(grid=xrgrid, frame=True)
fig.savefig("grdview.png")

fig = pygmt.Figure()
fig.grdimage(grid=xrgrid, frame=True)
fig.savefig("grdimage.png")

xrgrid.gmt.gtype is 1, so the grid is a geographic grid. However, the grid is still plotted as a Cartesian grid.

grdview grdimage
image image

System information

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

PyGMT information:
  version: v0.3.2.dev65+g9b57cc97.d20210404
System information:
  python: 3.9.2 | packaged by conda-forge | (default, Feb 21 2021, 05:02:20)  [Clang 11.0.1 ]
  executable: /Users/user/.miniconda/envs/pygmt/bin/python
  machine: macOS-11.2.3-x86_64-i386-64bit
Dependency information:
  numpy: 1.20.1
  pandas: 1.2.3
  xarray: 0.17.0
  netCDF4: 1.5.6
  packaging: 20.9
  ghostscript: 9.53.3
  gmt: 6.1.1
GMT library information:
  binary dir: /Users/user/.miniconda/envs/pygmt/bin
  cores: 8
  grid layout: rows
  library path: /Users/user/.miniconda/envs/pygmt/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/user/.miniconda/envs/pygmt/lib/gmt/plugins
  share dir: /Users/user/.miniconda/envs/pygmt/share/gmt
  version: 6.1.1
@seisman seisman added the bug Something isn't working label Apr 4, 2021
@seisman
Copy link
Member Author

seisman commented Apr 4, 2021

@PaulWessel This looks like an upstream GMT bug.

Here is a minimal script to reproduce the issue. It extracts a subregion from the earth_relief_01d grid, and then plot it using grdview. It's a geographic grid and should use a geographic projection, but it uses a Cartesian projection (see the image above).

import pygmt

xrgrid = pygmt.grdcut("@earth_relief_01d", region=[-116, -109, -47, -44])

fig = pygmt.Figure()
fig.grdview(grid=xrgrid, frame=True)
fig.savefig("grdview.png")

Internally, PyGMT calls the GMT_Create_Data function, and then calls GMT_Put_Matrix. Here are the arguments passed to GMT_Create_Data:

GMT_Create_Data(API, GMT_IS_GRID|GMT_VIA_MATRIX, GMT_IS_SURFACE, GMT_CONTAINER_ONLY|GMT_GRID_IS_GEO, ranges, incs, GMT_GRID_PIXEL_REG)

@PaulWessel
Copy link
Member

Here is what happens: Your GMT_GRID_IS_GEO bit-flag only leads to a call to gmt_set_geographic (GMT, GMT_OUT); inside GMT_Create_Data. Never mind that you did not pass GMT_IS_OUTPUT so this should be considered a grid for input, not output [I think I need to fix that bug which probably is not affecting many things anyway]. But even so, once your fig.grdview command is run a new session is created (right?) so there is no memory of setting anything to geographic. The solution is clearly to embed this info in the GMT_MATRIX structure, and that will mean in its void *hidden structure since there is no attribute in the GMT_MATRIX structure to hold that information. If we add a flag in the hidden struct that is set to geo then the GMT_Read_Data branch dealing with input matrices can explore that setting and call gmt_set_geographic as needed. I think that is the solution here. Agree?

@seisman
Copy link
Member Author

seisman commented Apr 5, 2021

The issue has been fixed in GenericMappingTools/gmt#5091 and GenericMappingTools/gmt#5092

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working upstream Bug or missing feature of upstream core GMT
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants