Wrap GMT's standard data type GMT_DATASET for table inputs

pygmt/clib/session.py

@@ -21,6 +21,7 @@
     vectors_to_arrays,
 )
 from pygmt.clib.loading import load_libgmt
+from pygmt.datatypes import GMT_DATASET
 from pygmt.exceptions import (
     GMTCLibError,
     GMTCLibNoSessionError,
@@ -1680,3 +1681,70 @@ def extract_region(self):
         if status != 0:
             raise GMTCLibError("Failed to extract region from current figure.")
         return wesn
+
+    def read_virtualfile(self, vfname):
+        """
+        Read data from a virtual file.
+
+        Parameters
+        ----------
+        vfname : str
+            Name of the virtual file to read.
+
+        Returns
+        -------
+        Pointer to the data, which can be casted into GMT data types.
+        """
+        c_read_virtualfile = self.get_libgmt_func(
+            "GMT_Read_VirtualFile",
+            argtypes=[
+                ctp.c_void_p,
+                ctp.c_char_p,
+            ],
+            restype=ctp.c_void_p,
+        )
+        return c_read_virtualfile(self.session_pointer, vfname.encode())
+
+    @contextmanager
+    def virtualfile_to_gmtdataset(self):
+        """
+        Create a virtual file for writing a GMT_DATASET object.
+
+        Yields
+        ------
+        vfile : str
+            Name of the virtual file.
+        """
+        family = "GMT_IS_DATASET"
+        geometry = "GMT_IS_PLP"
+        with self.open_virtual_file(family, geometry, "GMT_OUT", None) as vfile:
+            yield vfile
+
+    def gmtdataset_to_vectors(self, vfile):
+        """
+        Read GMT_DATASET object from a virtual file and convert to vectors.
+
+        Parameters
+        ----------
+        vfile : str
+            Name of the virtual file.
+
+        Returns
+        -------
+        vectors : list of 1-D arrays
+            List of vectors containing the data from the GMT_DATASET object.
+        """
+        # Read the virtual file and cast it to a pointer to a GMT_DATASET
+        ds = ctp.cast(self.read_virtualfile(vfile), ctp.POINTER(GMT_DATASET)).contents
+
+        # Loop over the tables, segments, and columns to get the data as vectors
+        vectors = []
+        for itbl in range(ds.n_tables):
+            dtbl = ds.table[itbl].contents
+            for iseg in range(dtbl.n_segments):
+                dseg = dtbl.segment[iseg].contents
+                for icol in range(dseg.n_columns):
+                    vectors.append(
+                        np.ctypeslib.as_array(dseg.data[icol], shape=(dseg.n_rows,))
+                    )
+        return vectors

pygmt/datatypes.py

@@ -0,0 +1,73 @@
+"""
+GMT data types for ctypes.
+
+See the GMT source code gmt_resources.h for the original C struct definitions.
+"""
+import ctypes as ctp
+
+
+class GMT_DATASEGMENT(ctp.Structure):
+    """
+    For holding segment lines in memory.
+    """
+
+    _fields_ = [
+        ("n_rows", ctp.c_uint64),  # Number of points in this segment
+        ("n_columns", ctp.c_uint64),  # Number of fields in each record (>= 2)
+        ("min", ctp.POINTER(ctp.c_double)),  # Minimum coordinate for each column
+        ("max", ctp.POINTER(ctp.c_double)),  # Maximum coordinate for each column
+        # Data x, y, and possibly other columns
+        ("data", ctp.POINTER(ctp.POINTER(ctp.c_double))),
+        ("label", ctp.c_char_p),  # Label string (if applicable)
+        ("header", ctp.c_char_p),  # Segment header (if applicable)
+        ("text", ctp.POINTER(ctp.c_char_p)),  # text beyond the data
+        ("hidden", ctp.c_void_p),  # Book-keeping variables "hidden" from the API
+    ]
+
+
+class GMT_DATATABLE(ctp.Structure):
+    """
+    To hold an array of line segment structures and header information in one
+    container.
+    """
+
+    _fields_ = [
+        ("n_headers", ctp.c_uint),  # Number of file header records (0 if no header)
+        ("n_columns", ctp.c_uint64),  # Number of columns (fields) in each record
+        ("n_segments", ctp.c_uint64),  # Number of segments in the array
+        ("n_records", ctp.c_uint64),  # Total number of data records across all segments
+        ("min", ctp.POINTER(ctp.c_double)),  # Minimum coordinate for each column
+        ("max", ctp.POINTER(ctp.c_double)),  # Maximum coordinate for each column
+        # Array with all file header records, if any
+        ("header", ctp.POINTER(ctp.c_char_p)),
+        # Pointer to array of segments
+        ("segment", ctp.POINTER(ctp.POINTER(GMT_DATASEGMENT))),
+        ("hidden", ctp.c_void_p),  # Book-keeping variables "hidden" from the API
+    ]
+
+
+class GMT_DATASET(ctp.Structure):
+    """
+    Single container for an array of GMT tables (files).
+    """
+
+    _fields_ = [
+        ("n_tables", ctp.c_uint64),  # The total number of tables (files) contained
+        ("n_columns", ctp.c_uint64),  # The number of data columns
+        ("n_segments", ctp.c_uint64),  #  The total number of segments across all tables
+        # The total number of data records across all tables
+        ("n_records", ctp.c_uint64),
+        ("min", ctp.POINTER(ctp.c_double)),  # Minimum coordinate for each column
+        ("max", ctp.POINTER(ctp.c_double)),  # Maximum coordinate for each column
+        # Pointer to array of tables
+        ("table", ctp.POINTER(ctp.POINTER(GMT_DATATABLE))),
+        # The datatype (numerical, text, or mixed) of this dataset
+        ("type", ctp.c_int32),
+        ("geometry", ctp.c_int32),  # The geometry of this dataset
+        # To store a referencing system string in PROJ.4 format
+        ("ProjRefPROJ4", ctp.c_char_p),
+        # To store a referencing system string in WKT format
+        ("ProjRefWKT", ctp.c_char_p),
+        ("ProjRefEPSG", ctp.c_int),  # To store a referencing system EPSG code
+        ("hidden", ctp.c_void_p),  # Book-keeping variables "hidden" from the API
+    ]

pygmt/src/grd2xyz.py

@@ -3,17 +3,12 @@
 """
 import warnings
 
+import numpy as np
 import pandas as pd
 import xarray as xr
 from pygmt.clib import Session
 from pygmt.exceptions import GMTInvalidInput
-from pygmt.helpers import (
-    GMTTempFile,
-    build_arg_string,
-    fmt_docstring,
-    kwargs_to_strings,
-    use_alias,
-)
+from pygmt.helpers import build_arg_string, fmt_docstring, kwargs_to_strings, use_alias
 
 __doctest_skip__ = ["grd2xyz"]
 
@@ -172,25 +167,37 @@
         # Reverse the dims because it is rows, columns ordered.
         dataframe_header = [grid.dims[1], grid.dims[0], grid.name]
 
-    with GMTTempFile() as tmpfile:
-        with Session() as lib:
-            file_context = lib.virtualfile_from_data(check_kind="raster", data=grid)
-            with file_context as infile:
-                if outfile is None:
-                    outfile = tmpfile.name
-                lib.call_module(
-                    module="grd2xyz",
-                    args=build_arg_string(kwargs, infile=infile, outfile=outfile),
-                )
-
-        # Read temporary csv output to a pandas table
-        if outfile == tmpfile.name:  # if user did not set outfile, return pd.DataFrame
-            result = pd.read_csv(
-                tmpfile.name, sep="\t", names=dataframe_header, comment=">"
+    with Session() as lib:
+        with lib.virtualfile_from_data(
+            check_kind="raster", data=grid
+        ) as invfile, lib.virtualfile_to_gmtdataset() as outvfile:
+            # Option 1
+            lib.call_module(
+                module="grd2xyz",
+                args=build_arg_string(kwargs, infile=invfile, outfile=outvfile),
+            )
+
+            if output_type == "file":
+                lib.call_module("write", f"{outvfile} {outfile} -Td")
+                return None
+            vectors = lib.gmtdataset_to_vectors(outvfile)
+            if output_type == "numpy":
+                return np.array(vectors).T
+            return pd.DataFrame(data=np.array(vectors).T, columns=dataframe_header)
+
+            """
+            # Option 2
+            if output_type == "file":
+                outvfile = outfile
+            lib.call_module(
+                module="grd2xyz",
+                args=build_arg_string(kwargs, infile=invfile, outfile=outvfile),
             )
-        elif outfile != tmpfile.name:  # return None if outfile set, output in outfile
-            result = None
 
-        if output_type == "numpy":
-            result = result.to_numpy()
-    return result
+            if output_type == "file":
+                return None
+            vectors = lib.gmtdataset_to_vectors(outvfile)
+            if output_type == "numpy":
+                return np.array(vectors).T
+            return pd.DataFrame(data=np.array(vectors).T, columns=dataframe_header)
+            """

pygmt/src/grdtrack.py

@@ -1,16 +1,11 @@
 """
 grdtrack - Sample grids at specified (x,y) locations.
 """
+import numpy as np
 import pandas as pd
 from pygmt.clib import Session
 from pygmt.exceptions import GMTInvalidInput
-from pygmt.helpers import (
-    GMTTempFile,
-    build_arg_string,
-    fmt_docstring,
-    kwargs_to_strings,
-    use_alias,
-)
+from pygmt.helpers import build_arg_string, fmt_docstring, kwargs_to_strings, use_alias
 
 __doctest_skip__ = ["grdtrack"]
 
@@ -43,7 +38,9 @@
     w="wrap",
 )
 @kwargs_to_strings(R="sequence", S="sequence", i="sequence_comma", o="sequence_comma")
-def grdtrack(grid, points=None, newcolname=None, outfile=None, **kwargs):
+def grdtrack(
+    grid, points=None, newcolname=None, output_type="pandas", outfile=None, **kwargs
+):
     r"""
     Sample grids at specified (x,y) locations.
 
@@ -292,29 +289,30 @@
     if hasattr(points, "columns") and newcolname is None:
         raise GMTInvalidInput("Please pass in a str to 'newcolname'")
 
-    with GMTTempFile(suffix=".csv") as tmpfile:
-        with Session() as lib:
-            with lib.virtualfile_from_data(
-                check_kind="raster", data=grid
-            ) as grdfile, lib.virtualfile_from_data(
-                check_kind="vector", data=points, required_data=False
-            ) as csvfile:
-                kwargs["G"] = grdfile
-                if outfile is None:  # Output to tmpfile if outfile is not set
-                    outfile = tmpfile.name
-                lib.call_module(
-                    module="grdtrack",
-                    args=build_arg_string(kwargs, infile=csvfile, outfile=outfile),
-                )
+    with Session() as lib:
+        with lib.virtualfile_from_data(
+            check_kind="raster", data=grid
+        ) as grdfile, lib.virtualfile_from_data(
+            check_kind="vector", data=points, required_data=False
+        ) as csvfile, lib.virtualfile_to_gmtdataset() as outvfile:
+            kwargs["G"] = grdfile
+            lib.call_module(
+                module="grdtrack",
+                args=build_arg_string(kwargs, infile=csvfile, outfile=outvfile),
+            )
+        if outfile is not None:
+            # if output_type == "file":
+            lib.call_module("write", f"{outvfile} {outfile} -Td")
+            return None
+
+        vectors = lib.gmtdataset_to_vectors(outvfile)
+
+        if output_type == "numpy":
+            return np.array(vectors).T
 
-        # Read temporary csv output to a pandas table
-        if outfile == tmpfile.name:  # if user did not set outfile, return pd.DataFrame
-            try:
-                column_names = points.columns.to_list() + [newcolname]
-                result = pd.read_csv(tmpfile.name, sep="\t", names=column_names)
-            except AttributeError:  # 'str' object has no attribute 'columns'
-                result = pd.read_csv(tmpfile.name, sep="\t", header=None, comment=">")
-        elif outfile != tmpfile.name:  # return None if outfile set, output in outfile
-            result = None
+        if isinstance(points, pd.DataFrame):
+            column_names = points.columns.to_list() + [newcolname]
+        else:
+            column_names = None
 
-    return result
+        return pd.DataFrame(np.array(vectors).T, columns=column_names)