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

[seisman]

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

[weiji14]

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.

Good spotting, this won't be specific to just grdcut, but the other grd* modulles as well. It will likely require some modification to these lines:

pygmt/pygmt/gridops.py

Lines 114 to 115 in f401f85

	with xr.open_dataarray(outgrid) as dataarray:
	result = dataarray.load()

There's 2 solutions I can see:

1. Trigger the GMTDataArrayAccessor inside the with block, so that the registration and gtype information is stored before the file disappears.
2. Don't delete the tempfile, maybe keep it somewhere while the session is still running.

Which would you prefer? Or is there a third way to resolve this?

[seisman]

Don't delete the tempfile, maybe keep it somewhere while the session is still running.

We probably don't want to maintain a list of temporary files.

Trigger the GMTDataArrayAccessor inside the with block, so that the registration and gtype information is stored before the file disappears.

It sounds a good way, but how to trigger it?

BTW, at least we need to check if the netCDF source exists in GMTDataArrayAccessor.

[weiji14]

It sounds a good way, but how to trigger it?

Simply by calling grid.gmt.registration inside the with block. That would then store the registration information in the class.

BTW, at least we need to check if the netCDF source exists in GMTDataArrayAccessor.

Yes, the error message could be improved, will see how best to do this in the PR.

[seisman]

Simply by calling grid.gmt.registration inside the with block. That would then store the registration information in the class.

Do you mean the following code?

 with xr.open_dataarray(outgrid) as dataarray: 
     result = dataarray.load() 
     result.gmt.registration
     result.gmt.gtype

It's a little weird.

[weiji14]

I think just result.gmt will do, but yes it's weird (and I'm sure pylint will complain about some useless-statement error). We'll also need to apply this to every grd* function that loads a temporary netcdf file. Maybe there's a better solution?

[seisman]

I think just result.gmt will do, but yes it's weird (and I'm sure pylint will complain about some useless-statement error). We'll also need to apply this to every grd* function that loads a temporary netcdf file. Maybe there's a better solution?

Perhaps we can add result.gmt as a temporary workaround and see how it can be improved in the future.

Issues #496 and #489 can't be fixed without the workaround.

[weiji14]

Ok, I'll try and do that tonight after dinner.

[seisman]

Keep the issue open, since we don't really like the workaround.

[seisman]

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

[weiji14]

Can't say I didn't see this one coming. I'm tempted to go with Solution 2 "Don't delete the tempfile, maybe keep it somewhere while the session is still running." now, since there's no easy way for to transfer the registration/gtype info from the original xarray grid to the sliced xarray grid.

Another option is to check explicitly if the *.nc file exists as you mentions. But if it doesn't, and we fallback on the default gridline registration/Cartesian setting, that's not ideal either. We'll need to come up with a better solution somehow.

[seisman]

The ultimate goal is to avoid temporary files. It's only possible if PyGMT supports GMT data types