|
| 1 | +""" |
| 2 | +Wrapper for the GMT_GRID_HEADER data structure and related utility functions. |
| 3 | +""" |
| 4 | + |
| 5 | +import ctypes as ctp |
| 6 | +from typing import Any, ClassVar |
| 7 | + |
| 8 | +import numpy as np |
| 9 | + |
| 10 | +# Constants for lengths of grid header variables. |
| 11 | +# |
| 12 | +# Note: Ideally we should be able to get these constants from the GMT shared library |
| 13 | +# using the ``lib["GMT_GRID_UNIT_LEN80"]`` syntax, but it causes cyclic import error. |
| 14 | +# So we have to hardcode the values here. |
| 15 | +GMT_GRID_UNIT_LEN80 = 80 |
| 16 | +GMT_GRID_TITLE_LEN80 = 80 |
| 17 | +GMT_GRID_VARNAME_LEN80 = 80 |
| 18 | +GMT_GRID_COMMAND_LEN320 = 320 |
| 19 | +GMT_GRID_REMARK_LEN160 = 160 |
| 20 | + |
| 21 | +# GMT uses single-precision for grids by default, but can be built to use |
| 22 | +# double-precision. Currently, only single-precision is supported. |
| 23 | +gmt_grdfloat = ctp.c_float |
| 24 | + |
| 25 | + |
| 26 | +def _parse_nameunits(nameunits: str) -> tuple[str, str | None]: |
| 27 | + """ |
| 28 | + Get the long_name and units attributes from x_units/y_units/z_units in the grid |
| 29 | + header. |
| 30 | +
|
| 31 | + In the GMT grid header, the x_units/y_units/z_units are strings in the form of |
| 32 | + ``long_name [units]``, in which both ``long_name`` and ``units`` are standard |
| 33 | + netCDF attributes defined by CF conventions. The ``[units]`` part is optional. |
| 34 | +
|
| 35 | + This function parses the x_units/y_units/z_units strings and gets the ``long_name`` |
| 36 | + and ``units`` attributes. |
| 37 | +
|
| 38 | + Parameters |
| 39 | + ---------- |
| 40 | + nameunits |
| 41 | + The x_units/y_units/z_units strings in the grid header. |
| 42 | +
|
| 43 | + Returns |
| 44 | + ------- |
| 45 | + (long_name, units) |
| 46 | + Tuple of netCDF attributes ``long_name`` and ``units``. ``units`` may be |
| 47 | + ``None``. |
| 48 | +
|
| 49 | + Examples |
| 50 | + -------- |
| 51 | + >>> _parse_nameunits("longitude [degrees_east]") |
| 52 | + ('longitude', 'degrees_east') |
| 53 | + >>> _parse_nameunits("latitude [degrees_north]") |
| 54 | + ('latitude', 'degrees_north') |
| 55 | + >>> _parse_nameunits("x") |
| 56 | + ('x', None) |
| 57 | + >>> _parse_nameunits("y") |
| 58 | + ('y', None) |
| 59 | + >>> |
| 60 | + """ |
| 61 | + parts = nameunits.split("[") |
| 62 | + long_name = parts[0].strip() |
| 63 | + units = parts[1].strip("]").strip() if len(parts) > 1 else None |
| 64 | + return long_name, units |
| 65 | + |
| 66 | + |
| 67 | +class _GMT_GRID_HEADER(ctp.Structure): # noqa: N801 |
| 68 | + """ |
| 69 | + GMT grid header structure for metadata about the grid. |
| 70 | +
|
| 71 | + The class is used in the `GMT_GRID`/`GMT_IMAGE`/`GMT_CUBE` data structure. See the |
| 72 | + GMT source code gmt_resources.h for the original C structure definitions. |
| 73 | + """ |
| 74 | + |
| 75 | + _fields_: ClassVar = [ |
| 76 | + # Number of columns |
| 77 | + ("n_columns", ctp.c_uint32), |
| 78 | + # Number of rows |
| 79 | + ("n_rows", ctp.c_uint32), |
| 80 | + # Grid registration, 0 for gridline and 1 for pixel |
| 81 | + ("registration", ctp.c_uint32), |
| 82 | + # Minimum/maximum x and y coordinates |
| 83 | + ("wesn", ctp.c_double * 4), |
| 84 | + # Minimum z value |
| 85 | + ("z_min", ctp.c_double), |
| 86 | + # Maximum z value |
| 87 | + ("z_max", ctp.c_double), |
| 88 | + # x and y increments |
| 89 | + ("inc", ctp.c_double * 2), |
| 90 | + # Grid values must be multiplied by this factor |
| 91 | + ("z_scale_factor", ctp.c_double), |
| 92 | + # After scaling, add this offset |
| 93 | + ("z_add_offset", ctp.c_double), |
| 94 | + # Units in x-directions, in the form "long_name [units]" |
| 95 | + ("x_units", ctp.c_char * GMT_GRID_UNIT_LEN80), |
| 96 | + # Units in y-direction, in the form "long_name [units]" |
| 97 | + ("y_units", ctp.c_char * GMT_GRID_UNIT_LEN80), |
| 98 | + # Grid value units, in the form "long_name [units]" |
| 99 | + ("z_units", ctp.c_char * GMT_GRID_UNIT_LEN80), |
| 100 | + # Name of data set |
| 101 | + ("title", ctp.c_char * GMT_GRID_TITLE_LEN80), |
| 102 | + # Name of generating command |
| 103 | + ("command", ctp.c_char * GMT_GRID_COMMAND_LEN320), |
| 104 | + # Comments for this data set |
| 105 | + ("remark", ctp.c_char * GMT_GRID_REMARK_LEN160), |
| 106 | + # Below are items used internally by GMT |
| 107 | + # Number of data points (n_columns * n_rows) [paddings are excluded] |
| 108 | + ("nm", ctp.c_size_t), |
| 109 | + # Actual number of items (not bytes) required to hold this grid (mx * my), |
| 110 | + # per band (for images) |
| 111 | + ("size", ctp.c_size_t), |
| 112 | + # Bits per data value (e.g., 32 for ints/floats; 8 for bytes). |
| 113 | + # Only used for ERSI ArcInfo ASCII Exchange grids. |
| 114 | + ("bits", ctp.c_uint), |
| 115 | + # For complex grid. |
| 116 | + # 0 for normal |
| 117 | + # GMT_GRID_IS_COMPLEX_REAL = real part of complex grid |
| 118 | + # GMT_GRID_IS_COMPLEX_IMAG = imag part of complex grid |
| 119 | + ("complex_mode", ctp.c_uint), |
| 120 | + # Grid format |
| 121 | + ("type", ctp.c_uint), |
| 122 | + # Number of bands [1]. Used with GMT_IMAGE containers |
| 123 | + ("n_bands", ctp.c_uint), |
| 124 | + # Actual x-dimension in memory. mx = n_columns + pad[0] + pad[1] |
| 125 | + ("mx", ctp.c_uint), |
| 126 | + # Actual y-dimension in memory. my = n_rows + pad[2] + pad[3] |
| 127 | + ("my", ctp.c_uint), |
| 128 | + # Paddings on west, east, south, north sides [2,2,2,2] |
| 129 | + ("pad", ctp.c_uint * 4), |
| 130 | + # Three or four char codes T|B R|C S|R|S (grd) or B|L|P + A|a (img) |
| 131 | + # describing array layout in mem and interleaving |
| 132 | + ("mem_layout", ctp.c_char * 4), |
| 133 | + # Missing value as stored in grid file |
| 134 | + ("nan_value", gmt_grdfloat), |
| 135 | + # 0.0 for gridline grids and 0.5 for pixel grids |
| 136 | + ("xy_off", ctp.c_double), |
| 137 | + # Referencing system string in PROJ.4 format |
| 138 | + ("ProjRefPROJ4", ctp.c_char_p), |
| 139 | + # Referencing system string in WKT format |
| 140 | + ("ProjRefWKT", ctp.c_char_p), |
| 141 | + # Referencing system EPSG code |
| 142 | + ("ProjRefEPSG", ctp.c_int), |
| 143 | + # Lower-level information for GMT use only |
| 144 | + ("hidden", ctp.c_void_p), |
| 145 | + ] |
| 146 | + |
| 147 | + def _parse_dimensions(self): |
| 148 | + """ |
| 149 | + Get dimension names and attributes from the grid header. |
| 150 | +
|
| 151 | + For a 2-D grid, the dimension names are set to "y" and "x" by default. The |
| 152 | + attributes for each dimension are parsed from the grid header following GMT |
| 153 | + source codes. See the GMT functions "gmtnc_put_units", "gmtnc_get_units" and |
| 154 | + "gmtnc_grd_info" for reference. |
| 155 | + """ |
| 156 | + # Default dimension names. |
| 157 | + dims = ("y", "x") |
| 158 | + nameunits = (self.y_units, self.x_units) |
| 159 | + |
| 160 | + # Dictionary for dimension attributes with the dimension name as the key. |
| 161 | + attrs = {dim: {} for dim in dims} |
| 162 | + # Dictionary for mapping the default dimension names to the actual names. |
| 163 | + newdims = {dim: dim for dim in dims} |
| 164 | + # Loop over dimensions and get the dimension name and attributes from header. |
| 165 | + for dim, nameunit in zip(dims, nameunits, strict=True): |
| 166 | + # The long_name and units attributes. |
| 167 | + long_name, units = _parse_nameunits(nameunit.decode()) |
| 168 | + if long_name: |
| 169 | + attrs[dim]["long_name"] = long_name |
| 170 | + if units: |
| 171 | + attrs[dim]["units"] = units |
| 172 | + |
| 173 | + # "degrees_east"/"degrees_north" are the units for geographic coordinates |
| 174 | + # following CF-conventions. |
| 175 | + if units == "degrees_east": |
| 176 | + attrs[dim]["standard_name"] = "longitude" |
| 177 | + newdims[dim] = "lon" |
| 178 | + elif units == "degrees_north": |
| 179 | + attrs[dim]["standard_name"] = "latitude" |
| 180 | + newdims[dim] = "lat" |
| 181 | + |
| 182 | + # Axis attributes are "X"/"Y"/"Z"/"T" for horizontal/vertical/time axis. |
| 183 | + attrs[dim]["axis"] = dim.upper() |
| 184 | + idx = 2 if dim == "y" else 0 |
| 185 | + attrs[dim]["actual_range"] = np.array(self.wesn[idx : idx + 2]) |
| 186 | + |
| 187 | + # Save the lists of dimension names and attributes in the _nc attribute. |
| 188 | + self._nc = { |
| 189 | + "dims": [newdims[dim] for dim in dims], |
| 190 | + "attrs": [attrs[dim] for dim in dims], |
| 191 | + } |
| 192 | + |
| 193 | + @property |
| 194 | + def name(self) -> str: |
| 195 | + """ |
| 196 | + Name of the grid. |
| 197 | + """ |
| 198 | + return "z" |
| 199 | + |
| 200 | + @property |
| 201 | + def data_attrs(self) -> dict[str, Any]: |
| 202 | + """ |
| 203 | + Attributes for the data variable from the grid header. |
| 204 | + """ |
| 205 | + attrs: dict[str, Any] = {} |
| 206 | + attrs["Conventions"] = "CF-1.7" |
| 207 | + attrs["title"] = self.title.decode() |
| 208 | + attrs["history"] = self.command.decode() |
| 209 | + attrs["description"] = self.remark.decode() |
| 210 | + long_name, units = _parse_nameunits(self.z_units.decode()) |
| 211 | + if long_name: |
| 212 | + attrs["long_name"] = long_name |
| 213 | + if units: |
| 214 | + attrs["units"] = units |
| 215 | + attrs["actual_range"] = np.array([self.z_min, self.z_max]) |
| 216 | + return attrs |
| 217 | + |
| 218 | + @property |
| 219 | + def dims(self) -> list: |
| 220 | + """ |
| 221 | + List of dimension names. |
| 222 | + """ |
| 223 | + if not hasattr(self, "_nc"): |
| 224 | + self._parse_dimensions() |
| 225 | + return self._nc["dims"] |
| 226 | + |
| 227 | + @property |
| 228 | + def dim_attrs(self) -> list[dict]: |
| 229 | + """ |
| 230 | + List of attributes for each dimension. |
| 231 | + """ |
| 232 | + if not hasattr(self, "_nc"): |
| 233 | + self._parse_dimensions() |
| 234 | + return self._nc["attrs"] |
| 235 | + |
| 236 | + @property |
| 237 | + def gtype(self) -> int: |
| 238 | + """ |
| 239 | + Grid type. 0 for Cartesian grid and 1 for geographic grid. |
| 240 | +
|
| 241 | + The grid is assumed to be Cartesian by default. If the x/y dimensions are named |
| 242 | + "lon"/"lat" or have units "degrees_east"/"degrees_north", then the grid is |
| 243 | + assumed to be geographic. |
| 244 | + """ |
| 245 | + dims = self.dims |
| 246 | + gtype = 1 if dims[0] == "lat" and dims[1] == "lon" else 0 |
| 247 | + return gtype |
0 commit comments