Use read_virtualfile_to_data to hide the technical details

Author: seisman

pygmt/clib/session.py

@@ -10,6 +10,7 @@
 import warnings
 from contextlib import contextmanager, nullcontext
 
+from pygmt.datatypes import GMT_GRID
 import numpy as np
 import pandas as pd
 from packaging.version import Version
@@ -1647,6 +1648,28 @@ def read_virtualfile(self, vfname):
         )
         return c_read_virtualfile(self.session_pointer, vfname.encode())
 
+    def read_virtualfile_to_data(self, vfname, kind):
+        """
+        Read a virtual file and convert to a GMT data container.
+
+        Parameters
+        ----------
+        vfname : str
+            Name of the virtual file to read.
+        kind : str
+            The kind of data container to create. Choose from "grid" or
+            "dataset".
+
+        Returns
+        -------
+        Pointer to the GMT_GRID or GMT_DATASET data container.
+        """
+        type = {
+            "grid": GMT_GRID,
+            # "dataset": GMT_DATASET, # implemented in PR #2729
+        }[kind]
+        return ctp.cast(self.read_virtualfile(vfname), ctp.POINTER(type))
+
     @contextmanager
     def virtualfile_to_data(self, kind):
         """
@@ -1703,8 +1726,6 @@ def _add_pad(x):
                 np.zeros((2, x.shape[1] + 4)),
             ]
 
-        from pygmt.datatypes import GMT_GRID
-
         family = "GMT_IS_GRID"
         geometry = "GMT_IS_SURFACE"
         matrix, region, inc = dataarray_to_matrix(xrgrid)

pygmt/datatypes.py

@@ -7,14 +7,22 @@
 
 import numpy as np
 import xarray as xr
-from pygmt.clib.session import Session
-
-# Lengths of grid header variables
-with Session() as lib:
-    GMT_GRID_UNIT_LEN80 = lib["GMT_GRID_UNIT_LEN80"]
-    GMT_GRID_TITLE_LEN80 = lib["GMT_GRID_TITLE_LEN80"]
-    GMT_GRID_COMMAND_LEN320 = lib["GMT_GRID_COMMAND_LEN320"]
-    GMT_GRID_REMARK_LEN160 = lib["GMT_GRID_REMARK_LEN160"]
+
+# from pygmt.clib.session import Session
+
+# # Lengths of grid header variables
+# with Session() as lib:
+#     GMT_GRID_UNIT_LEN80 = lib["GMT_GRID_UNIT_LEN80"]
+#     GMT_GRID_TITLE_LEN80 = lib["GMT_GRID_TITLE_LEN80"]
+#     GMT_GRID_COMMAND_LEN320 = lib["GMT_GRID_COMMAND_LEN320"]
+#     GMT_GRID_REMARK_LEN160 = lib["GMT_GRID_REMARK_LEN160"]
+# Ideally we should be able to get the lengths of grid header variables
+# from GMT shared library, but it cause cyclic import error.
+# So we hardcode the values here.
+GMT_GRID_UNIT_LEN80 = 80
+GMT_GRID_TITLE_LEN80 = 80
+GMT_GRID_COMMAND_LEN320 = 320
+GMT_GRID_REMARK_LEN160 = 160
 
 # GMT uses single-precision for grids by default, but can be built to use
 # double-precision. Currently, only single-precision is supported.

pygmt/src/grdcut.py

@@ -109,5 +109,6 @@ def grdcut(grid, outgrid=None, **kwargs):
             if outgrid is not None:
                 lib.call_module("write", f"{outfile} {outgrid} -Tg")
                 return None
-            gmtgrid = lib.read_virtualfile(outfile)
-            return ctp.cast(gmtgrid, ctp.POINTER(GMT_GRID)).contents.to_dataarray()
+
+            gmtgrid = lib.read_virtualfile_to_data(outfile, kind="grid")
+            return gmtgrid.contents.to_dataarray()