From 048ccc929145a703a48e5dd0aed6acd624b76010 Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Thu, 5 Apr 2018 17:43:20 -1000 Subject: [PATCH 1/5] Create LibGMT.grid_to_vfile for xarray support Implemented a doctest using the earth relief data. Not passing at the moment. Failing with a Seg Fault. Might be because the direction needs to be `GMT_IN|GMT_IS_REFERENCE`. Need support for | in the direction argument. --- gmt/clib/core.py | 123 +++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 108 insertions(+), 15 deletions(-) diff --git a/gmt/clib/core.py b/gmt/clib/core.py index 0516fb345a5..032d26e773f 100644 --- a/gmt/clib/core.py +++ b/gmt/clib/core.py @@ -845,6 +845,7 @@ def vectors_to_vfile(self, *vectors): Examples -------- + >>> from gmt.helpers import GMTTempFile >>> import numpy as np >>> import pandas as pd >>> x = [1, 2, 3] @@ -853,14 +854,11 @@ def vectors_to_vfile(self, *vectors): >>> with LibGMT() as lib: ... with lib.vectors_to_vfile(x, y, z) as vfile: ... # Send the output to a file so that we can read it - ... ofile = 'vectors_to_vfile_example.txt' - ... lib.call_module('info', '{} ->{}'.format(vfile, ofile)) - >>> with open(ofile) as f: - ... # Replace tabs with spaces - ... print(f.read().strip().replace('\\t', ' ')) + ... with GMTTempFile() as ofile: + ... args = '{} ->{}'.format(vfile, ofile.name) + ... lib.call_module('info', args) + ... print(ofile.read().strip()) : N = 3 <1/3> <4/6> <7/9> - >>> # Clean up the output file - >>> os.remove(ofile) """ arrays = vectors_to_arrays(vectors) @@ -894,7 +892,6 @@ def matrix_to_vfile(self, matrix): The virtual file will contain the array as a GMT Dataset (via matrix). - **Not meant for creating GMT Grids**. The grid requires more metadata than just the data matrix. This creates a Dataset (table). @@ -925,6 +922,7 @@ def matrix_to_vfile(self, matrix): Examples -------- + >>> from gmt.helpers import GMTTempFile >>> import numpy as np >>> data = np.arange(12).reshape((4, 3)) >>> print(data) @@ -935,14 +933,11 @@ def matrix_to_vfile(self, matrix): >>> with LibGMT() as lib: ... with lib.matrix_to_vfile(data) as vfile: ... # Send the output to a file so that we can read it - ... ofile = 'matrix_to_vfile_example.txt' - ... lib.call_module('info', '{} ->{}'.format(vfile, ofile)) - >>> with open(ofile) as f: - ... # Replace tabs with spaces - ... print(f.read().strip().replace('\\t', ' ')) + ... with GMTTempFile() as ofile: + ... args = '{} ->{}'.format(vfile, ofile.name) + ... lib.call_module('info', args) + ... print(ofile.read().strip()) : N = 4 <0/9> <1/10> <2/11> - >>> # Clean up the output file - >>> os.remove(ofile) """ rows, columns = matrix.shape @@ -963,6 +958,104 @@ def matrix_to_vfile(self, matrix): with self.open_virtual_file(*vf_args) as vfile: yield vfile + @contextmanager + def grid_to_vfile(self, grid): + """ + Store a grid in a GMT virtual file with a GMT_GRID container. + + Used to pass grid data into GMT modules. Grids must be + ``xarray.DataArray`` instances. + + Context manager (use in a ``with`` block). Yields the virtual file name + that you can pass as an argument to a GMT module call. Closes the + virtual file upon exit of the ``with`` block. + + The virtual file will contain the grid as a ``GMT_GRID`` (via matrix). + + Use this instead of creating ``GMT_GRID`` and virtual files by hand + with :meth:`~gmt.clib.LibGMT.create_data`, + :meth:`~gmt.clib.LibGMT.put_matrix`, and + :meth:`~gmt.clib.LibGMT.open_virtual_file` + + The grid data matrix must be C contiguous in memory. If it is not + (e.g., it is a slice of a larger array), the array will be copied to + make sure it is. + + Parameters + ---------- + grid : xarray.DataArraw + The grid that will be included in the virtual file. + + Yields + ------ + vfile : str + The name of virtual file. Pass this as a file name argument to a + GMT module. + + Examples + -------- + + >>> import xarray as xr + >>> from gmt.datasets import load_earth_relief + >>> from gmt.helpers import GMTTempFile + >>> data = load_earth_relief(resolution='60m') + >>> print(data.shape) + (181, 361) + >>> print(data.lon.values.min(), data.lon.values.max()) + -180.0 180.0 + >>> print(data.lat.values.min(), data.lat.values.max()) + -90.0 90.0 + >>> print(data.values.min(), data.values.max()) + -8425.0 5551.0 + >>> with LibGMT() as lib: + ... with lib.grid_to_vfile(data) as vfile: + ... # Send the output to a file so that we can read it + ... with GMTTempFile() as ofile: + ... args = '{} -L0 ->{}'.format(vfile, ofile.name) + ... lib.call_module('grdinfo', args) + ... print(ofile.read().strip()) + : N = 4 <0/9> <1/10> <2/11> + + """ + if len(grid.dims) != 2: + raise GMTInvalidInput( + "Invalid number of grid dimensions '{}'. Must be 2." + .format(len(grid.dims))) + + # Get the region and grid increment from grid coordinates. + # dims is ordered as row, column. + y, x = [grid.coords[dim].values for dim in grid.dims] + region = [x.min(), x.max(), y.min(), y.max()] + x_incs = x[1:] - x[0:-1] + x_inc = x_incs[0] + if not np.allclose(x_incs, x_inc): + raise GMTInvalidInput( + "Grid appears to have irregular spacing in the '{}' dimension." + .format(grid.dims[1])) + y_incs = y[1:] - y[0:-1] + y_inc = y_incs[0] + + nrows, ncolumns = grid.shape + + matrix = grid.values[:] + # Array must be in C contiguous order to pass its memory pointer to GMT + if not matrix.flags.c_contiguous: + matrix = matrix.copy(order='C') + + family = 'GMT_IS_GRID|GMT_VIA_MATRIX' + geometry = 'GMT_IS_SURFACE' + + gmt_grid = self.create_data(family, geometry, + mode='GMT_CONTAINER_ONLY', + ranges=region, inc=[x_inc, y_inc], + dim=[ncolumns, nrows, 1, 0]) + + self.put_matrix(gmt_grid, matrix) + + args = (family, geometry, 'GMT_IN|GMT_IS_REFERENCE', gmt_grid) + with self.open_virtual_file(*args) as vfile: + yield vfile + def extract_region(self): """ Extract the WESN bounding box of the currently active figure. From 1c5fb67f8ad495ea038b883b80c55e0486e7af9c Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Fri, 6 Apr 2018 17:39:25 -1000 Subject: [PATCH 2/5] Breakup the function into more testable parts Still crashing when using ranges and inc but works with dim. Might be a bug in GMT --- gmt/clib/core.py | 53 +++++++++++++++-------------------------------- gmt/clib/utils.py | 46 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 36 deletions(-) diff --git a/gmt/clib/core.py b/gmt/clib/core.py index 032d26e773f..0a6a92baf33 100644 --- a/gmt/clib/core.py +++ b/gmt/clib/core.py @@ -9,7 +9,8 @@ import numpy as np from ..exceptions import GMTCLibError, GMTCLibNoSessionError, GMTInvalidInput -from .utils import load_libgmt, kwargs_to_ctypes_array, vectors_to_arrays +from .utils import load_libgmt, kwargs_to_ctypes_array, vectors_to_arrays, \ + dataarray_to_matrix class LibGMT(): # pylint: disable=too-many-instance-attributes @@ -408,9 +409,7 @@ def create_data(self, family, geometry, mode, **kwargs): # Parse and check input arguments family_int = self._parse_constant(family, valid=self.data_families, valid_modifiers=self.data_vias) - if mode not in self.data_modes: - raise GMTInvalidInput( - "Invalid data creation mode '{}'.".format(mode)) + mode_int = self._parse_constant(mode, valid=self.data_modes) geometry_int = self._parse_constant(geometry, valid=self.data_geometries) # dim is required (get a segmentation fault if passing it as None @@ -430,7 +429,7 @@ def create_data(self, family, geometry, mode, **kwargs): # container should be created empty. Fill it in later using put_vector # and put_matrix. data_ptr = c_create_data(self.current_session, family_int, - geometry_int, self.get_constant(mode), dim, + geometry_int, mode_int, dim, ranges, inc, registration, pad, None) if data_ptr is None: @@ -534,6 +533,11 @@ def _check_dtype_and_dim(self, array, ndim): ... gmttype = lib._check_dtype_and_dim(data, ndim=1) ... gmttype == lib.get_constant('GMT_DOUBLE') True + >>> data = np.ones((5, 2), dtype='float32') + >>> with LibGMT() as lib: + ... gmttype = lib._check_dtype_and_dim(data, ndim=2) + ... gmttype == lib.get_constant('GMT_FLOAT') + True """ if array.dtype.name not in self._dtypes: @@ -789,9 +793,10 @@ def open_virtual_file(self, family, geometry, direction, data): valid_modifiers=self.data_vias) geometry_int = self._parse_constant(geometry, valid=self.data_geometries) - if direction not in ('GMT_IN', 'GMT_OUT'): - raise GMTInvalidInput("Invalid direction '{}'.".format(direction)) - direction_int = self.get_constant(direction) + direction_int = self._parse_constant( + direction, valid=['GMT_IN', 'GMT_OUT'], + valid_modifiers=['GMT_IS_REFERENCE']) + # valid_modifiers=None) buff = ctypes.create_string_buffer(self.get_constant('GMT_STR16')) @@ -995,7 +1000,6 @@ def grid_to_vfile(self, grid): Examples -------- - >>> import xarray as xr >>> from gmt.datasets import load_earth_relief >>> from gmt.helpers import GMTTempFile >>> data = load_earth_relief(resolution='60m') @@ -1014,46 +1018,23 @@ def grid_to_vfile(self, grid): ... args = '{} -L0 ->{}'.format(vfile, ofile.name) ... lib.call_module('grdinfo', args) ... print(ofile.read().strip()) - : N = 4 <0/9> <1/10> <2/11> """ - if len(grid.dims) != 2: - raise GMTInvalidInput( - "Invalid number of grid dimensions '{}'. Must be 2." - .format(len(grid.dims))) - - # Get the region and grid increment from grid coordinates. - # dims is ordered as row, column. - y, x = [grid.coords[dim].values for dim in grid.dims] - region = [x.min(), x.max(), y.min(), y.max()] - x_incs = x[1:] - x[0:-1] - x_inc = x_incs[0] - if not np.allclose(x_incs, x_inc): - raise GMTInvalidInput( - "Grid appears to have irregular spacing in the '{}' dimension." - .format(grid.dims[1])) - y_incs = y[1:] - y[0:-1] - y_inc = y_incs[0] - + region, inc, matrix = dataarray_to_matrix(grid) nrows, ncolumns = grid.shape - - matrix = grid.values[:] - # Array must be in C contiguous order to pass its memory pointer to GMT - if not matrix.flags.c_contiguous: - matrix = matrix.copy(order='C') - family = 'GMT_IS_GRID|GMT_VIA_MATRIX' geometry = 'GMT_IS_SURFACE' gmt_grid = self.create_data(family, geometry, mode='GMT_CONTAINER_ONLY', - ranges=region, inc=[x_inc, y_inc], - dim=[ncolumns, nrows, 1, 0]) + dim=[ncolumns, nrows, 1, 0], + ranges=region, inc=inc) self.put_matrix(gmt_grid, matrix) args = (family, geometry, 'GMT_IN|GMT_IS_REFERENCE', gmt_grid) with self.open_virtual_file(*args) as vfile: + # yield 'bla' yield vfile def extract_region(self): diff --git a/gmt/clib/utils.py b/gmt/clib/utils.py index 7561aaa0b04..aec8769deda 100644 --- a/gmt/clib/utils.py +++ b/gmt/clib/utils.py @@ -10,6 +10,52 @@ from ..exceptions import GMTOSError, GMTCLibError, GMTCLibNotFoundError +def dataarray_to_matrix(grid): + """ + Transform a xarray.DataArray into region, increment, and data matrix. + + Examples + -------- + + >>> from gmt.datasets import load_earth_relief + >>> grid = load_earth_relief(resolution='60m') + >>> region, inc, matrix = grid_to_region_inc_matrix(grid) + >>> print(region) + [-180.0, 180.0, -90.0, 90.0] + >>> print(inc) + [1.0, 1.0] + >>> type(matrix) + + >>> print(matrix.shape) + (181, 361) + + """ + if len(grid.dims) != 2: + raise GMTInvalidInput( + "Invalid number of grid dimensions '{}'. Must be 2." + .format(len(grid.dims))) + + # Get the region and grid increment from grid coordinates. + # dims is ordered as row, column. + y, x = [grid.coords[dim].values for dim in grid.dims] + region = [x.min(), x.max(), y.min(), y.max()] + x_incs = x[1:] - x[0:-1] + x_inc = x_incs[0] + if not np.allclose(x_incs, x_inc): + raise GMTInvalidInput( + "Grid appears to have irregular spacing in the '{}' dimension." + .format(grid.dims[1])) + y_incs = y[1:] - y[0:-1] + y_inc = y_incs[0] + + matrix = grid.values + # Array must be in C contiguous order to pass its memory pointer to GMT + if not matrix.flags.c_contiguous: + matrix = matrix.copy(order='C') + + return region, [x_inc, y_inc], matrix + + def vectors_to_arrays(vectors): """ Convert 1d vectors (lists, arrays or pandas.Series) to C contiguous 1d From 4a5a85f76bbe88f67e05435e7e4a60152a473ec0 Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Mon, 9 Apr 2018 10:40:33 -1000 Subject: [PATCH 3/5] Tests pass on GMT trunk Need docs for dataarray_to_matrix and more tests for the code. --- gmt/clib/core.py | 12 ++++-------- gmt/clib/utils.py | 23 +++++++++++------------ 2 files changed, 15 insertions(+), 20 deletions(-) diff --git a/gmt/clib/core.py b/gmt/clib/core.py index 0a6a92baf33..d5212ca05f2 100644 --- a/gmt/clib/core.py +++ b/gmt/clib/core.py @@ -1015,26 +1015,22 @@ def grid_to_vfile(self, grid): ... with lib.grid_to_vfile(data) as vfile: ... # Send the output to a file so that we can read it ... with GMTTempFile() as ofile: - ... args = '{} -L0 ->{}'.format(vfile, ofile.name) + ... args = '{} -L0 -Cn ->{}'.format(vfile, ofile.name) ... lib.call_module('grdinfo', args) ... print(ofile.read().strip()) + -180 180 -90 90 -8425 5551 1 1 361 181 + >>> # The output is: w e s n z0 z1 dx dy n_columns n_rows """ - region, inc, matrix = dataarray_to_matrix(grid) - nrows, ncolumns = grid.shape + matrix, region, inc = dataarray_to_matrix(grid) family = 'GMT_IS_GRID|GMT_VIA_MATRIX' geometry = 'GMT_IS_SURFACE' - gmt_grid = self.create_data(family, geometry, mode='GMT_CONTAINER_ONLY', - dim=[ncolumns, nrows, 1, 0], ranges=region, inc=inc) - self.put_matrix(gmt_grid, matrix) - args = (family, geometry, 'GMT_IN|GMT_IS_REFERENCE', gmt_grid) with self.open_virtual_file(*args) as vfile: - # yield 'bla' yield vfile def extract_region(self): diff --git a/gmt/clib/utils.py b/gmt/clib/utils.py index aec8769deda..61aa27cb0f0 100644 --- a/gmt/clib/utils.py +++ b/gmt/clib/utils.py @@ -19,7 +19,7 @@ def dataarray_to_matrix(grid): >>> from gmt.datasets import load_earth_relief >>> grid = load_earth_relief(resolution='60m') - >>> region, inc, matrix = grid_to_region_inc_matrix(grid) + >>> matrix, region, inc = dataarray_to_matrix(grid) >>> print(region) [-180.0, 180.0, -90.0, 90.0] >>> print(inc) @@ -34,26 +34,25 @@ def dataarray_to_matrix(grid): raise GMTInvalidInput( "Invalid number of grid dimensions '{}'. Must be 2." .format(len(grid.dims))) - # Get the region and grid increment from grid coordinates. - # dims is ordered as row, column. - y, x = [grid.coords[dim].values for dim in grid.dims] - region = [x.min(), x.max(), y.min(), y.max()] - x_incs = x[1:] - x[0:-1] - x_inc = x_incs[0] - if not np.allclose(x_incs, x_inc): + east, north = [grid.coords[dim].values for dim in grid.dims] + region = [north.min(), north.max(), east.min(), east.max()] + north_incs = north[1:] - north[0:-1] + north_inc = north_incs[0] + if not np.allclose(north_incs, north_inc): raise GMTInvalidInput( "Grid appears to have irregular spacing in the '{}' dimension." .format(grid.dims[1])) - y_incs = y[1:] - y[0:-1] - y_inc = y_incs[0] + east_incs = east[1:] - east[0:-1] + east_inc = east_incs[0] + inc = [north_inc, east_inc] - matrix = grid.values + matrix = grid.values[:] # Array must be in C contiguous order to pass its memory pointer to GMT if not matrix.flags.c_contiguous: matrix = matrix.copy(order='C') - return region, [x_inc, y_inc], matrix + return matrix, region, inc def vectors_to_arrays(vectors): From 555a8ce0522b409e8395027fc87069e66132cc7c Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Mon, 9 Apr 2018 18:15:50 -1000 Subject: [PATCH 4/5] Implement tests and docs for all new functions --- gmt/clib/core.py | 87 ++++++++++++++++++++++------------- gmt/clib/utils.py | 101 +++++++++++++++++++++++++++++++---------- gmt/tests/test_clib.py | 32 +++++++++++-- 3 files changed, 160 insertions(+), 60 deletions(-) diff --git a/gmt/clib/core.py b/gmt/clib/core.py index d5212ca05f2..e366d9d0e9f 100644 --- a/gmt/clib/core.py +++ b/gmt/clib/core.py @@ -10,7 +10,7 @@ from ..exceptions import GMTCLibError, GMTCLibNoSessionError, GMTInvalidInput from .utils import load_libgmt, kwargs_to_ctypes_array, vectors_to_arrays, \ - dataarray_to_matrix + dataarray_to_matrix, as_c_contiguous class LibGMT(): # pylint: disable=too-many-instance-attributes @@ -78,6 +78,11 @@ class LibGMT(): # pylint: disable=too-many-instance-attributes 'GMT_OUTPUT', ] + grid_registrations = [ + 'GMT_GRID_PIXEL_REG', + 'GMT_GRID_NODE_REG', + ] + # Map numpy dtypes to GMT types _dtypes = { 'float64': 'GMT_DOUBLE', @@ -378,9 +383,9 @@ def create_data(self, family, geometry, mode, **kwargs): The increments between points of the dataset. See the C function documentation. registration : int - The node registration (what the coordinates mean). Also a very - complex argument. See the C function documentation. Defaults to - ``GMT_GRID_NODE_REG``. + The node registration (what the coordinates mean). Can be + ``'GMT_GRID_PIXEL_REG'`` or ``'GMT_GRID_NODE_REG'``. Defaults to + ``'GMT_GRID_NODE_REG'``. pad : int The grid padding. Defaults to ``GMT_PAD_DEFAULT``. @@ -406,31 +411,35 @@ def create_data(self, family, geometry, mode, **kwargs): ] c_create_data.restype = ctypes.c_void_p - # Parse and check input arguments family_int = self._parse_constant(family, valid=self.data_families, valid_modifiers=self.data_vias) mode_int = self._parse_constant(mode, valid=self.data_modes) - geometry_int = self._parse_constant(geometry, - valid=self.data_geometries) - # dim is required (get a segmentation fault if passing it as None - if 'dim' not in kwargs: - kwargs['dim'] = [0]*4 - # Convert dim, ranges, and inc to ctypes arrays if given + geometry_int = self._parse_constant( + geometry, valid=self.data_geometries) + registration_int = self._parse_constant( + kwargs.get('registration', 'GMT_GRID_NODE_REG'), + valid=self.grid_registrations) + + # Convert dim, ranges, and inc to ctypes arrays if given (will be None + # if not given to represent NULL pointers) dim = kwargs_to_ctypes_array('dim', kwargs, ctypes.c_uint64*4) ranges = kwargs_to_ctypes_array('ranges', kwargs, ctypes.c_double*4) inc = kwargs_to_ctypes_array('inc', kwargs, ctypes.c_double*2) - pad = self._parse_pad(family, kwargs) - - # Use the GMT defaults if no value is given - registration = kwargs.get('registration', - self.get_constant('GMT_GRID_NODE_REG')) # Use a NULL pointer (None) for existing data to indicate that the # container should be created empty. Fill it in later using put_vector # and put_matrix. - data_ptr = c_create_data(self.current_session, family_int, - geometry_int, mode_int, dim, - ranges, inc, registration, pad, None) + data_ptr = c_create_data( + self.current_session, + family_int, + geometry_int, + mode_int, + dim, + ranges, + inc, + registration_int, + self._parse_pad(family, kwargs), + None) if data_ptr is None: raise GMTCLibError("Failed to create an empty GMT data pointer.") @@ -750,6 +759,7 @@ def open_virtual_file(self, family, geometry, direction, data): Examples -------- + >>> from gmt.helpers import GMTTempFile >>> import os >>> import numpy as np >>> x = np.array([0, 1, 2, 3, 4]) @@ -768,15 +778,12 @@ def open_virtual_file(self, family, geometry, direction, data): ... # Add the dataset to a virtual file ... vfargs = (family, geometry, 'GMT_IN', dataset) ... with lib.open_virtual_file(*vfargs) as vfile: - ... # Send the output to a file so that we can read it - ... ofile = 'virtual_file_example.txt' - ... lib.call_module('info', '{} ->{}'.format(vfile, ofile)) - >>> with open(ofile) as f: - ... # Replace tabs with spaces - ... print(f.read().strip().replace('\\t', ' ')) + ... # Send the output to a temp file so that we can read it + ... with GMTTempFile() as ofile: + ... args = '{} ->{}'.format(vfile, ofile.name) + ... lib.call_module('info', args) + ... print(ofile.read().strip()) : N = 5 <0/4> <5/9> - >>> # Clean up the output file - >>> os.remove(ofile) """ c_open_virtualfile = self._libgmt.GMT_Open_VirtualFile @@ -796,7 +803,6 @@ def open_virtual_file(self, family, geometry, direction, data): direction_int = self._parse_constant( direction, valid=['GMT_IN', 'GMT_OUT'], valid_modifiers=['GMT_IS_REFERENCE']) - # valid_modifiers=None) buff = ctypes.create_string_buffer(self.get_constant('GMT_STR16')) @@ -866,6 +872,13 @@ def vectors_to_vfile(self, *vectors): : N = 3 <1/3> <4/6> <7/9> """ + # Conversion to a C-contiguous array needs to be done here and not in + # put_matrix because we need to maintain a reference to the copy while + # it is being used by the C API. Otherwise, the array would be garbage + # collected and the memory freed. Creating it in this context manager + # guarantees that the copy will be around until the virtual file is + # closed. + # The conversion is implicit in vectors_to_arrays. arrays = vectors_to_arrays(vectors) columns = len(arrays) @@ -945,11 +958,14 @@ def matrix_to_vfile(self, matrix): : N = 4 <0/9> <1/10> <2/11> """ + # Conversion to a C-contiguous array needs to be done here and not in + # put_matrix because we need to maintain a reference to the copy while + # it is being used by the C API. Otherwise, the array would be garbage + # collected and the memory freed. Creating it in this context manager + # guarantees that the copy will be around until the virtual file is + # closed. + matrix = as_c_contiguous(matrix) rows, columns = matrix.shape - # Array must be in C contiguous order to pass its memory pointer to - # GMT. - if not matrix.flags.c_contiguous: - matrix = matrix.copy(order='C') family = 'GMT_IS_DATASET|GMT_VIA_MATRIX' geometry = 'GMT_IS_POINT' @@ -1022,6 +1038,13 @@ def grid_to_vfile(self, grid): >>> # The output is: w e s n z0 z1 dx dy n_columns n_rows """ + # Conversion to a C-contiguous array needs to be done here and not in + # put_matrix because we need to maintain a reference to the copy while + # it is being used by the C API. Otherwise, the array would be garbage + # collected and the memory freed. Creating it in this context manager + # guarantees that the copy will be around until the virtual file is + # closed. + # The conversion is implicit in dataarray_to_matrix. matrix, region, inc = dataarray_to_matrix(grid) family = 'GMT_IS_GRID|GMT_VIA_MATRIX' geometry = 'GMT_IS_SURFACE' diff --git a/gmt/clib/utils.py b/gmt/clib/utils.py index 61aa27cb0f0..cb555ebbe22 100644 --- a/gmt/clib/utils.py +++ b/gmt/clib/utils.py @@ -7,17 +7,48 @@ import numpy as np import pandas -from ..exceptions import GMTOSError, GMTCLibError, GMTCLibNotFoundError +from ..exceptions import GMTOSError, GMTCLibError, GMTCLibNotFoundError, \ + GMTInvalidInput def dataarray_to_matrix(grid): """ - Transform a xarray.DataArray into region, increment, and data matrix. + Transform a xarray.DataArray into a data 2D array and metadata. + + Use this to extract the underlying numpy array of data and the region and + increment for the grid. + + Only allows grids with two dimensions and constant grid spacing (GMT + doesn't allow variable grid spacing). + + If the underlying data array is not C contiguous, for example if it's a + slice of a larger grid, a copy will need to be generated. + + Parameters + ---------- + grid : xarray.DataArray + The input grid as a DataArray instance. Information is retrieved from + the coordinate arrays, not from headers. + + Returns + ------- + matrix : 2d-array + The 2D array of data from the grid. + region : list + The West, East, South, North boundaries of the grid. + inc : list + The grid spacing in East-West and North-South, respectively. + + Raises + ------ + GMTInvalidInput + If the grid has more than two dimensions or variable grid spacing. Examples -------- >>> from gmt.datasets import load_earth_relief + >>> # Use the global Earth relief grid with 1 degree spacing (60') >>> grid = load_earth_relief(resolution='60m') >>> matrix, region, inc = dataarray_to_matrix(grid) >>> print(region) @@ -28,30 +59,52 @@ def dataarray_to_matrix(grid): >>> print(matrix.shape) (181, 361) + >>> matrix.flags.c_contiguous + True + >>> # Using a slice of the grid, the matrix will be copied to guarantee + >>> # that it's C-contiguous in memory. The increment should be unchanged. + >>> matrix, region, inc = dataarray_to_matrix(grid[10:41,30:101]) + >>> matrix.flags.c_contiguous + True + >>> print(matrix.shape) + (31, 71) + >>> print(region) + [-150.0, -80.0, -80.0, -50.0] + >>> print(inc) + [1.0, 1.0] + >>> # but not if only taking every other grid point. + >>> matrix, region, inc = dataarray_to_matrix(grid[10:41:2,30:101:2]) + >>> matrix.flags.c_contiguous + True + >>> print(matrix.shape) + (16, 36) + >>> print(region) + [-150.0, -80.0, -80.0, -50.0] + >>> print(inc) + [2.0, 2.0] """ if len(grid.dims) != 2: raise GMTInvalidInput( "Invalid number of grid dimensions '{}'. Must be 2." .format(len(grid.dims))) - # Get the region and grid increment from grid coordinates. - east, north = [grid.coords[dim].values for dim in grid.dims] - region = [north.min(), north.max(), east.min(), east.max()] - north_incs = north[1:] - north[0:-1] - north_inc = north_incs[0] - if not np.allclose(north_incs, north_inc): - raise GMTInvalidInput( - "Grid appears to have irregular spacing in the '{}' dimension." - .format(grid.dims[1])) - east_incs = east[1:] - east[0:-1] - east_inc = east_incs[0] - inc = [north_inc, east_inc] - - matrix = grid.values[:] - # Array must be in C contiguous order to pass its memory pointer to GMT - if not matrix.flags.c_contiguous: - matrix = matrix.copy(order='C') - + # Extract region and inc from the grid + region = [] + inc = [] + # Reverse the dims because it is rows, columns ordered. In geographic + # grids, this would be North-South, East-West. GMT's region and inc are + # East-West, North-South. + for dim in grid.dims[::-1]: + coord = grid.coords[dim].values + coord_incs = coord[1:] - coord[0:-1] + coord_inc = coord_incs[0] + if not np.allclose(coord_incs, coord_inc): + raise GMTInvalidInput( + "Grid appears to have irregular spacing in the '{}' dimension." + .format(dim)) + region.extend([coord.min(), coord.max()]) + inc.append(coord_inc) + matrix = as_c_contiguous(grid.values[:]) return matrix, region, inc @@ -98,11 +151,11 @@ def vectors_to_arrays(vectors): True """ - arrays = [_as_c_contiguous(_as_array(i)) for i in vectors] + arrays = [as_c_contiguous(_as_array(i)) for i in vectors] return arrays -def _as_c_contiguous(array): +def as_c_contiguous(array): """ Ensure a numpy array is C contiguous in memory. @@ -128,7 +181,7 @@ def _as_c_contiguous(array): array([1, 3, 5]) >>> x.flags.c_contiguous False - >>> new_x = _as_c_contiguous(x) + >>> new_x = as_c_contiguous(x) >>> new_x array([1, 3, 5]) >>> new_x.flags.c_contiguous @@ -136,7 +189,7 @@ def _as_c_contiguous(array): >>> x = np.array([8, 9, 10]) >>> x.flags.c_contiguous True - >>> _as_c_contiguous(x).flags.c_contiguous + >>> as_c_contiguous(x).flags.c_contiguous True """ diff --git a/gmt/tests/test_clib.py b/gmt/tests/test_clib.py index c847dd8fdfb..19dedfbac83 100644 --- a/gmt/tests/test_clib.py +++ b/gmt/tests/test_clib.py @@ -8,10 +8,12 @@ import pytest import numpy as np import numpy.testing as npt -import pandas +import pandas as pd +import xarray as xr from ..clib.core import LibGMT -from ..clib.utils import clib_extension, load_libgmt, check_libgmt +from ..clib.utils import clib_extension, load_libgmt, check_libgmt, \ + dataarray_to_matrix from ..exceptions import GMTCLibError, GMTOSError, GMTCLibNotFoundError, \ GMTCLibNoSessionError, GMTInvalidInput from ..helpers import GMTTempFile @@ -439,7 +441,7 @@ def test_put_matrix_grid(): mode='GMT_CONTAINER_ONLY', ranges=wesn[:4], inc=inc, - registration=lib.get_constant('GMT_GRID_NODE_REG'), + registration='GMT_GRID_NODE_REG', ) data = np.arange(shape[0]*shape[1], dtype=dtype).reshape(shape) lib.put_matrix(grid, matrix=data) @@ -610,7 +612,7 @@ def test_vectors_to_vfile_pandas(): dtypes = 'float32 float64 int32 int64 uint32 uint64'.split() size = 13 for dtype in dtypes: - data = pandas.DataFrame( + data = pd.DataFrame( data=dict(x=np.arange(size, dtype=dtype), y=np.arange(size, size*2, 1, dtype=dtype), z=np.arange(size*2, size*3, 1, dtype=dtype)) @@ -690,3 +692,25 @@ def test_write_data_fails(): with pytest.raises(GMTCLibError): lib.write_data('GMT_IS_VECTOR', 'GMT_IS_POINT', 'GMT_WRITE_SET', [1]*6, 'some-file-name', None) + + +def test_dataarray_to_matrix_dims_fails(): + "Check that it fails for > 2 dims" + # Make a 3D regular grid + data = np.ones((10, 12, 11), dtype='float32') + x = np.arange(11) + y = np.arange(12) + z = np.arange(10) + grid = xr.DataArray(data, coords=[('z', z), ('y', y), ('x', x)]) + with pytest.raises(GMTInvalidInput): + dataarray_to_matrix(grid) + + +def test_dataarray_to_matrix_inc_fails(): + "Check that it fails for variable increments" + data = np.ones((4, 5), dtype='float64') + x = np.linspace(0, 1, 5) + y = np.logspace(2, 3, 4) + grid = xr.DataArray(data, coords=[('y', y), ('x', x)]) + with pytest.raises(GMTInvalidInput): + dataarray_to_matrix(grid) From b3e8c8883ab7ec14bde565f22a9bc3b29ff63e8f Mon Sep 17 00:00:00 2001 From: Leonardo Uieda Date: Mon, 9 Apr 2018 18:22:16 -1000 Subject: [PATCH 5/5] Better docs for the _to_vfile methods --- gmt/clib/core.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/gmt/clib/core.py b/gmt/clib/core.py index e366d9d0e9f..10836bc323a 100644 --- a/gmt/clib/core.py +++ b/gmt/clib/core.py @@ -825,7 +825,7 @@ def open_virtual_file(self, family, geometry, direction, data): @contextmanager def vectors_to_vfile(self, *vectors): """ - Store 1d arrays in a GMT virtual file to pass them to a module. + Store 1d arrays in a GMT virtual file to use as a module input. Context manager (use in a ``with`` block). Yields the virtual file name that you can pass as an argument to a GMT module call. Closes the @@ -836,10 +836,11 @@ def vectors_to_vfile(self, *vectors): :meth:`~gmt.clib.LibGMT.put_vector`, and :meth:`~gmt.clib.LibGMT.open_virtual_file` - The virtual file will contain the arrays as a GMT Dataset (via - vectors). If the arrays are C contiguous blocks of memory, they will be - passed without copying to GMT. If they are not (e.g., they are columns - of a 2D array), they will need to be copied to a contiguous block. + The virtual file will contain the arrays as ``GMT Vector`` structures. + + If the arrays are C contiguous blocks of memory, they will be passed + without copying to GMT. If they are not (e.g., they are columns of a 2D + array), they will need to be copied to a contiguous block. Parameters ---------- @@ -902,13 +903,13 @@ def vectors_to_vfile(self, *vectors): @contextmanager def matrix_to_vfile(self, matrix): """ - Store a 2d array in a GMT virtual file Dataset to pass it to a module. + Store a 2d array in a GMT virtual file to use as a module input. Context manager (use in a ``with`` block). Yields the virtual file name that you can pass as an argument to a GMT module call. Closes the virtual file upon exit of the ``with`` block. - The virtual file will contain the array as a GMT Dataset (via matrix). + The virtual file will contain the array as a ``GMT_MATRIX``. **Not meant for creating GMT Grids**. The grid requires more metadata than just the data matrix. This creates a Dataset (table). @@ -982,7 +983,7 @@ def matrix_to_vfile(self, matrix): @contextmanager def grid_to_vfile(self, grid): """ - Store a grid in a GMT virtual file with a GMT_GRID container. + Store a grid in a GMT virtual file to use as a module input. Used to pass grid data into GMT modules. Grids must be ``xarray.DataArray`` instances. @@ -991,7 +992,7 @@ def grid_to_vfile(self, grid): that you can pass as an argument to a GMT module call. Closes the virtual file upon exit of the ``with`` block. - The virtual file will contain the grid as a ``GMT_GRID`` (via matrix). + The virtual file will contain the grid as a ``GMT_MATRIX``. Use this instead of creating ``GMT_GRID`` and virtual files by hand with :meth:`~gmt.clib.LibGMT.create_data`,