pygmt: GMTDataArrayAccessor doesn't work for temporary files that have been sliced

Description of the problem

It may be a GMTDataArrayAccessor bug, or maybe we should update grdcut() to work with GMTDataArrayAccessor. The problem is that the temporary netCDF file is already deleted when GMTDataArrayAccessor tries to access it.

Full code that generated the error

import pygmt

grid = pygmt.grdcut("@earth_relief_01d", region=[30, 50, 60, 70])
print(grid)
print(grid.gmt.registration)

Full error message

<xarray.DataArray 'z' (lat: 10, lon: 20)>
array([[  46. ,   -2.5,    0. ,   42. ,  150. ,  220. ,  172. ,  147. ,
         153.5,  159. ,  177. ,  172. ,  165. ,  137. ,  144.5,  150. ,
         119. ,  142. ,  131.5,  149.5],
       [  18.5,   19. ,   80.5,  139. ,  120.5,   33. ,   54. ,  202.5,
         153. ,  126. ,  175. ,  165.5,   69. ,  180. ,  147. ,  147. ,
          61. ,  106. ,  112.5,  147. ],
       [ 125.5,  173. ,  185. ,  136.5,   70.5,   51.5,  158.5,  140. ,
         142. ,   97. ,  139. ,  133. ,   58.5,   53.5,  109. ,  171. ,
         158. ,  153. ,  152. ,  156. ],
       [ 173.5,  223.5,  199.5,  130.5,   95. ,  130.5,  170.5,  114. ,
          51. ,   80. ,   78.5,   30.5,   82. ,  129. ,  137. ,  126. ,
         139. ,  118. ,  138.5,  148. ],
       [ 216.5,  163. ,  133. ,  121. ,   51.5,   -7. ,  -18. ,   65. ,
          96. ,   34.5,   19.5,   58.5,   78. ,   42. ,   69.5,   96.5,
          82. ,  155. ,  132.5,  192.5],
       [ 187. ,  163.5,  123.5,   74.5,   12. ,  -50. , -177.5, -126. ,
        -106.5,  -37. ,   94. ,  101.5,  121.5,   28. ,   38. ,   64.5,
          67. ,   89. ,  118.5,  154. ],
       [ 169. ,  137. ,   63.5,  -23. , -100. ,   33. ,   85. ,  151. ,
         179.5,  184.5,   60. ,  -34. ,    1. ,   -8.5,   10. ,   40.5,
          31. ,   33. ,   76.5,  136.5],
       [ 288. ,  252. ,  190. ,  273. ,  208.5,  231. ,  191.5,  207.5,
         233.5,  243. ,  180.5,  -38. ,  -27. ,  -10.5,   16.5,  -18. ,
         -28. ,  -31. ,   89. ,   59.5],
       [ 157.5,  138.5,  234. ,  218.5,  260.5,  214. ,  235. ,  171. ,
         -49.5,  -83. ,  -74. ,  -68. ,  -68. ,  -25. ,   46.5,   16.5,
         -43.5,  -51. ,  -31. ,  -26.5],
       [ 204.5,  145.5,   77.5, -152. , -122. , -178. , -187.5, -160.5,
        -125.5, -144.5, -154. , -138. ,  -90.5,  -65.5,  -81. ,  -76.5,
         -64.5,  -46.5,  -12. ,   -0.5]], dtype=float32)
Coordinates:
  * lon      (lon) float64 30.5 31.5 32.5 33.5 34.5 ... 45.5 46.5 47.5 48.5 49.5
  * lat      (lat) float64 60.5 61.5 62.5 63.5 64.5 65.5 66.5 67.5 68.5 69.5
Attributes:
    long_name:     elevation (m)
    actual_range:  [-187.5  288. ]
grdinfo [ERROR]: File /var/folders/mz/d3wtwjtx4fg8qlds5lgt1c240000gn/T/pygmt-6f1lhxvo.nc was not found
grdinfo [ERROR]: Cannot find file /var/folders/mz/d3wtwjtx4fg8qlds5lgt1c240000gn/T/pygmt-6f1lhxvo.nc
grdinfo [ERROR]: Must specify one or more input files
Traceback (most recent call last):
  File "test.py", line 5, in <module>
    print(grid.gmt.registration)
  File "/Users/seisman/.anaconda/lib/python3.7/site-packages/xarray/core/extensions.py", line 36, in __get__
    accessor_obj = self._accessor(obj)
  File "/Users/seisman/Gits/gmt/pygmt/pygmt/modules.py", line 247, in __init__
    int, grdinfo(self._source, C="n", o="10,11").split()
  File "/Users/seisman/Gits/gmt/pygmt/pygmt/modules.py", line 51, in grdinfo
    lib.call_module("grdinfo", arg_str)
  File "/Users/seisman/Gits/gmt/pygmt/pygmt/clib/session.py", line 502, in call_module
    module, status, self._error_message
pygmt.exceptions.GMTCLibError: Module 'grdinfo' failed with status code 71:
grdinfo [ERROR]: File /var/folders/mz/d3wtwjtx4fg8qlds5lgt1c240000gn/T/pygmt-6f1lhxvo.nc was not found
grdinfo [ERROR]: Cannot find file /var/folders/mz/d3wtwjtx4fg8qlds5lgt1c240000gn/T/pygmt-6f1lhxvo.nc
grdinfo [ERROR]: Must specify one or more input files

System information

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

PyGMT information:
  version: v0.1.2+11.g168c77c
System information:
  python: 3.7.6 (default, Jan  8 2020, 13:42:34)  [Clang 4.0.1 (tags/RELEASE_401/final)]
  executable: /Users/seisman/.anaconda/bin/python
  machine: Darwin-19.5.0-x86_64-i386-64bit
Dependency information:
  numpy: 1.18.1
  pandas: 1.0.1
  xarray: 0.15.1
  netCDF4: 1.5.3
  packaging: 20.1
  ghostscript: 9.52
  gmt: 6.2.0_308386e_2020.07.11
GMT library information:
  binary dir: /Users/seisman/.anaconda/bin
  cores: 4
  grid layout: rows
  library path: /Users/seisman/local/GMT/lib/libgmt.dylib
  padding: 2
  plugin dir: /Users/seisman/local/GMT/lib/gmt/plugins
  share dir: /Users/seisman/local/GMT/share
  version: 6.2.0

About this issue

Original URL
State: open
Created 4 years ago
Comments: 16 (16 by maintainers)

Most upvoted comments

For xarray.DataArray objects, currently PyGMT relies on the output of the grdinfo call to determine the registration and gtype. If the grdinfo fails, PyGMT fallbacks to the registration=0 and gtype=0.

See what GMT does at https://github.com/GenericMappingTools/gmt/blob/e63b2590b68969908efc393d008a9d024617eed2/src/gmt_nc.c#L656.

It seems GMT has its own logic to detect a netCDF file’s grid registration. I think PyGMT should follow the same logic. It would be easier if GMT can provide an API function for this.

seisman on Jul 16, 2022

It still doesn’t work for a slice of a grid. For example:

In [1]: import pygmt

In [2]: grid = pygmt.grdcut("@earth_relief_01d", region="0/90/0/90")

In [3]: grid.gmt.registration
Out[3]: 1

In [4]: grid2 = grid[10:41,30:101]

In [5]: grid2
Out[5]:
<xarray.DataArray 'z' (lat: 31, lon: 60)>
array([[  482. ,   435. ,   391.5, ..., -3402.5, -3340.5, -3281. ],
       [  731. ,   571. ,   393. , ..., -3303.5, -3228.5, -3176. ],
       [  543. ,   490. ,   395.5, ..., -3166.5, -3132. , -3066.5],
       ...,
       [ 1312.5,  1136.5,  1058. , ...,  1514.5,  3121. ,  3689.5],
       [ 1104.5,   993.5,  1084.5, ...,   875.5,   815.5,   894.5],
       [  623. ,  1156. ,  1329. , ...,   861.5,   836. ,   822. ]],
      dtype=float32)
Coordinates:
  * lon      (lon) float64 30.5 31.5 32.5 33.5 34.5 ... 85.5 86.5 87.5 88.5 89.5
  * lat      (lat) float64 10.5 11.5 12.5 13.5 14.5 ... 36.5 37.5 38.5 39.5 40.5
Attributes:
    long_name:     elevation (m)
    actual_range:  [-5122.   5651.5]

In [6]: grid2.gmt.registration
grdinfo [ERROR]: File /var/folders/mz/d3wtwjtx4fg8qlds5lgt1c240000gn/T/pygmt-15alnpr7.nc was not found
grdinfo [ERROR]: Cannot find file /var/folders/mz/d3wtwjtx4fg8qlds5lgt1c240000gn/T/pygmt-15alnpr7.nc
grdinfo [ERROR]: Must specify one or more input files

seisman on Jul 22, 2020