|
| 1 | +""" |
| 2 | +Wrapper for the GMT_CUBE data type. |
| 3 | +""" |
| 4 | + |
| 5 | +import ctypes as ctp |
| 6 | +from typing import ClassVar |
| 7 | + |
| 8 | +import numpy as np |
| 9 | +import xarray as xr |
| 10 | +from pygmt.datatypes.header import ( |
| 11 | + _GMT_GRID_HEADER, |
| 12 | + GMT_GRID_UNIT_LEN80, |
| 13 | + GMT_GRID_VARNAME_LEN80, |
| 14 | + _parse_nameunits, |
| 15 | + gmt_grdfloat, |
| 16 | +) |
| 17 | + |
| 18 | + |
| 19 | +class _GMT_CUBE(ctp.Structure): # noqa: N801 |
| 20 | + """ |
| 21 | + GMT cube data structure for 3D data. |
| 22 | + """ |
| 23 | + |
| 24 | + _fields_: ClassVar = [ |
| 25 | + # Pointer to full GMT 2-D header for a layer (common to all layers) |
| 26 | + ("header", ctp.POINTER(_GMT_GRID_HEADER)), |
| 27 | + # Pointer to the gmt_grdfloat 3-D cube - a stack of 2-D padded grids |
| 28 | + ("data", ctp.POINTER(gmt_grdfloat)), |
| 29 | + # Vector of x coordinates common to all layers |
| 30 | + ("x", ctp.POINTER(ctp.c_double)), |
| 31 | + # Vector of y coordinates common to all layers |
| 32 | + ("y", ctp.POINTER(ctp.c_double)), |
| 33 | + # Low-level information for GMT use only |
| 34 | + ("hidden", ctp.c_void_p), |
| 35 | + # GMT_CUBE_IS_STACK if input dataset was a list of 2-D grids rather than a |
| 36 | + # single cube |
| 37 | + ("mode", ctp.c_uint), |
| 38 | + # Minimum/max z values (complements header->wesn[4]) |
| 39 | + ("z_range", ctp.c_double * 2), |
| 40 | + # z increment (complements inc[2]) (0 if variable z spacing) |
| 41 | + ("z_inc", ctp.c_double), |
| 42 | + # Array of z values (complements x, y) |
| 43 | + ("z", ctp.POINTER(ctp.c_double)), |
| 44 | + # Name of the 3-D variable, if read from file (or empty if just one) |
| 45 | + ("name", ctp.c_char * GMT_GRID_VARNAME_LEN80), |
| 46 | + # Units in 3rd direction (complements x_units, y_units, z_units) |
| 47 | + ("units", ctp.c_char * GMT_GRID_UNIT_LEN80), |
| 48 | + ] |
| 49 | + |
| 50 | + def to_dataarray(self): |
| 51 | + """ |
| 52 | + Convert the GMT_CUBE to an xarray.DataArray. |
| 53 | +
|
| 54 | + Returns |
| 55 | + ------- |
| 56 | + xarray.DataArray: The data array representation of the GMT_CUBE. |
| 57 | + """ |
| 58 | + # The grid header |
| 59 | + header = self.header.contents |
| 60 | + |
| 61 | + name = "cube" |
| 62 | + # Dimensions and attributes |
| 63 | + dims = header.dims |
| 64 | + dim_attrs = header.dim_attrs |
| 65 | + |
| 66 | + # Patch for the 3rd dimension |
| 67 | + dims.append("z") |
| 68 | + z_attrs = {"actual_range": np.array(self.z_range[:]), "axis": "Z"} |
| 69 | + long_name, units = _parse_nameunits(self.units.decode()) |
| 70 | + if long_name: |
| 71 | + z_attrs["long_name"] = long_name |
| 72 | + if units: |
| 73 | + z_attrs["units"] = units |
| 74 | + dim_attrs.append(z_attrs) |
| 75 | + |
| 76 | + # The coordinates, given as a tuple of the form (dims, data, attrs) |
| 77 | + coords = [ |
| 78 | + (dims[0], self.y[: header.n_rows], dim_attrs[0]), |
| 79 | + (dims[1], self.x[: header.n_columns], dim_attrs[1]), |
| 80 | + # header->n_bands is used for the number of layers for 3-D cubes |
| 81 | + (dims[2], self.z[: header.n_bands], dim_attrs[1]), |
| 82 | + ] |
| 83 | + |
| 84 | + # The data array without paddings |
| 85 | + pad = header.pad[:] |
| 86 | + data = np.reshape( |
| 87 | + self.data[: header.mx * header.my * header.n_bands], |
| 88 | + (header.my, header.mx, header.n_bands), |
| 89 | + )[pad[2] : header.my - pad[3], pad[0] : header.mx - pad[1], :] |
| 90 | + |
| 91 | + # Create the xarray.DataArray object |
| 92 | + grid = xr.DataArray(data, coords=coords, name=name, attrs=header.dataA_attrs) |
| 93 | + return grid |
0 commit comments