Source code for harmonica._io.oasis_montaj_grd

# Copyright (c) 2018 The Harmonica Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
# The support for compressed grids was inspired by the code in
# https://github.com/Loop3D/geosoft_grid copyrighted by Loop3D and released
# under the MIT License.
"""
Function to read Oasis Montaj© .grd file
"""

import array
import zlib

import numpy as np
import xarray as xr

# Define the valid element sizes (ES variable) for GRD files
# (values > 1024 correspond to compressed versions of the grid)
VALID_ELEMENT_SIZES = (1, 2, 4, 8, 1024 + 1, 1024 + 2, 1024 + 4, 1024 + 8)
# Define dummy values for each data type
DUMMIES = {
    "b": -127,
    "B": 255,
    "h": -32767,
    "H": 65535,
    "i": -2147483647,
    "I": 4294967295,
    "f": -1e32,
    "d": -1e32,
}


[docs]def load_oasis_montaj_grid(fname): """ Reads gridded data from an Oasis Montaj© .grd file. The version 2 of the Geosoft© Grid File Format (GRD) stores gridded products in binary data. This function can read those files and parse the information in the header. It returns the data in a :class:`xarray.DataArray` for convenience. .. warning:: This function has not been tested against a wide range of GRD files. This could lead to incorrect readings of the stored data. Please report any unwanted behaviour by opening an issue in Harmonica: https://github.com/fatiando/harmonica/issues .. important:: This function only supports reading GRD files using the **version 2** of the Geosoft© Grid File Format. .. important:: This function is not supporting orderings different than ±1, or colour grids. Parameters ---------- fname : string or file-like object Path to the .grd file. Returns ------- grid : :class:`xarray.DataArray` :class:`xarray.DataArray` containing the grid, its coordinates and header information. References ---------- https://help.seequent.com/Oasis-montaj/9.9/en/Content/ss/glossary/grid_file_format__grd.htm https://github.com/Loop3D/geosoft_grid """ # Read the header and the grid array with open(fname, "rb") as grd_file: # Read the header (first 512 bytes) header = _read_header(grd_file.read(512)) # Check for valid flags _check_ordering(header["ordering"]) _check_sign_flag(header["sign_flag"]) # Get data type for the grid elements data_type = _get_data_type(header["n_bytes_per_element"], header["sign_flag"]) # Read grid grid = grd_file.read() # Decompress grid if needed if header["n_bytes_per_element"] > 1024: grid = _decompress_grid(grid) # Load the grid values as an array with the proper data_type grid = array.array(data_type, grid) # Convert to numpy array as float64 grid = np.array(grid, dtype=np.float64) # Remove dummy values grid = _remove_dummies(grid, data_type) # Scale the grid grid = np.array(grid / header["data_factor"] + header["base_value"]) # Reshape the grid based on the ordering if header["ordering"] == 1: order = "C" shape = (header["shape_v"], header["shape_e"]) spacing = (header["spacing_v"], header["spacing_e"]) elif header["ordering"] == -1: order = "F" shape = (header["shape_e"], header["shape_v"]) spacing = (header["spacing_e"], header["spacing_v"]) grid = grid.reshape(shape, order=order) # Build coords if header["rotation"] == 0: easting, northing = _build_coordinates( header["x_origin"], header["y_origin"], shape, spacing ) dims = ("northing", "easting") coords = {"easting": easting, "northing": northing} else: easting, northing = _build_rotated_coordinates( header["x_origin"], header["y_origin"], shape, spacing, header["rotation"] ) dims = ("y", "x") coords = {"easting": (dims, easting), "northing": (dims, northing)} # Build an xarray.DataArray for the grid grid = xr.DataArray( grid, coords=coords, dims=dims, attrs=header, ) return grid
def _read_header(header_bytes): """ Read GRD file header Parameters ---------- header_bytes : byte A sequence of 512 bytes containing the header of a GRD file. Returns ------- header : dict Dictionary containing the information present in the header. Notes ----- The GRD header consists in 512 contiguous bytes. It's divided in four sections: * Data Storage * Geographic Information * Data (Z) Scaling * Undefined Application Parameters """ header = {} # Read data storage ES, SF, NE, NV, KX = array.array("i", header_bytes[0 : 5 * 4]) # noqa: N806 header.update( { "n_bytes_per_element": ES, "sign_flag": SF, "shape_e": NE, "shape_v": NV, "ordering": KX, } ) # Read geographic info DE, DV, X0, Y0, ROT = array.array("d", header_bytes[20 : 20 + 5 * 8]) # noqa: N806 header.update( { "spacing_e": DE, "spacing_v": DV, "x_origin": X0, "y_origin": Y0, "rotation": ROT, } ) # Read data scaling ZBASE, ZMULT = array.array("d", header_bytes[60 : 60 + 2 * 8]) # noqa: N806 header.update( { "base_value": ZBASE, "data_factor": ZMULT, } ) # Read optional parameters # (ignore map LABEL and MAPNO) PROJ, UNITX, UNITY, UNITZ, NVPTS = array.array( # noqa: N806 "i", header_bytes[140 : 140 + 5 * 4] ) IZMIN, IZMAX, IZMED, IZMEA = array.array( # noqa: N806 "f", header_bytes[160 : 160 + 4 * 4] ) (ZVAR,) = array.array("d", header_bytes[176 : 176 + 8]) # noqa: N806 (PRCS,) = array.array("i", header_bytes[184 : 184 + 4]) # noqa: N806 header.update( { "map_projection": PROJ, "units_x": UNITX, "units_y": UNITY, "units_z": UNITZ, "n_valid_points": NVPTS, "grid_min": IZMIN, "grid_max": IZMAX, "grid_median": IZMED, "grid_mean": IZMEA, "grid_variance": ZVAR, "process_flag": PRCS, } ) return header def _check_ordering(ordering): """ Check if the ordering value is within the ones we are supporting """ if ordering not in (-1, 1): raise NotImplementedError( f"Found an ordering (a.k.a as KX) equal to '{ordering}'. " + "Only orderings equal to 1 and -1 are supported." ) def _check_sign_flag(sign_flag): """ Check if sign_flag value is within the ones we are supporting """ if sign_flag == 3: raise NotImplementedError( "Reading .grd files with colour grids is not currenty supported." ) def _get_data_type(n_bytes_per_element, sign_flag): """ Return the data type for the grid values References ---------- https://docs.python.org/3/library/array.html """ # Check if number of bytes per element is valid if n_bytes_per_element not in VALID_ELEMENT_SIZES: raise NotImplementedError( "Found a 'Grid data element size' (a.k.a. 'ES') value " + f"of '{n_bytes_per_element}'. " "Only values equal to 1, 2, 4 and 8 are valid, " + "along with their compressed counterparts (1025, 1026, 1028, 1032)." ) # Shift the n_bytes_per_element in case of compressed grids if n_bytes_per_element > 1024: n_bytes_per_element -= 1024 # Determine the data type of the grid elements if n_bytes_per_element == 1: if sign_flag == 0: data_type = "B" # unsigned char elif sign_flag == 1: data_type = "b" # signed char elif n_bytes_per_element == 2: if sign_flag == 0: data_type = "H" # unsigned short elif sign_flag == 1: data_type = "h" # signed short elif n_bytes_per_element == 4: if sign_flag == 0: data_type = "I" # unsigned int elif sign_flag == 1: data_type = "i" # signed int elif sign_flag == 2: data_type = "f" # float elif n_bytes_per_element == 8: data_type = "d" return data_type def _remove_dummies(grid, data_type): """ Replace dummy values for NaNs """ # Create dictionary with dummy value for each data type if data_type in ("f", "d"): grid[grid <= DUMMIES[data_type]] = np.nan return grid grid[grid == DUMMIES[data_type]] = np.nan return grid def _decompress_grid(grid_compressed): """ Decompress the grid using gzip Even if the header specifies that the grid is compressed using a LZRW1 algorithm, it's using gzip instead. The first two 4 bytes sequences correspond to the compression signature and to the compression type. We are going to ignore those and start reading from the number of blocks (offset 8). Parameters ---------- grid_compressed : bytes Sequence of bytes corresponding to the compressed grid. They should be every byte starting from offset 512 of the GRD file until its end. Returns ------- grid : bytes Uncompressed version of the ``grid_compressed`` parameter. """ # Number of blocks (n_blocks,) = array.array("i", grid_compressed[8 : 8 + 4]) # Number of vectors per block (vectors_per_block,) = array.array("i", grid_compressed[12 : 12 + 4]) # File offset from start of every block block_offsets = array.array("q", grid_compressed[16 : 16 + n_blocks * 8]) # Compressed size of every block compressed_block_sizes = array.array( "i", grid_compressed[16 + n_blocks * 8 : 16 + n_blocks * 8 + n_blocks * 4], ) # Combine grid grid = b"" # Read each block for i in range(n_blocks): # Define the start and end offsets for each compressed blocks # We need to remove the 512 to account for the missing header. # There is an unexplained 16 byte header that we also need to remove. start_offset = block_offsets[i] - 512 + 16 end_offset = compressed_block_sizes[i] + block_offsets[i] - 512 # Decompress the block grid_sub = zlib.decompress( grid_compressed[start_offset:end_offset], bufsize=zlib.DEF_BUF_SIZE, ) # Add it to the running grid grid += grid_sub return grid def _build_coordinates(west, south, shape, spacing): """ Create the coordinates for the grid Generates 1d arrays for the easting and northing coordinates of the grid. Assumes unrotated grids. Parameters ---------- west : float Westernmost coordinate of the grid. south : float Southernmost coordinate of the grid. shape : tuple Tuple of ints containing the number of elements along each direction in the following order: ``n_northing``, ``n_easting`` spacing : tuple Tuple of floats containing the distance between adjacent grid elements along each direction in the following order: ``s_northing``, ``s_easting``. Returns ------- easting : 1d-array Array containing the values of the easting coordinates of the grid. northing : 1d-array Array containing the values of the northing coordinates of the grid. """ easting = np.linspace(west, west + spacing[1] * (shape[1] - 1), shape[1]) northing = np.linspace(south, south + spacing[0] * (shape[0] - 1), shape[0]) return easting, northing def _build_rotated_coordinates(west, south, shape, spacing, rotation_deg): """ Create the coordinates for a rotated grid Generates 2d arrays for the easting and northing coordinates of the grid. Parameters ---------- west : float Westernmost coordinate of the grid. south : float Southernmost coordinate of the grid. shape : tuple Tuple of ints containing the number of elements along each unrotated direction in the following order: ``n_y``, ``n_x`` spacing : tuple Tuple of floats containing the distance between adjacent grid elements along each unrotated direction in the following order: ``spacing_y``, ``spacing_x``. rotation_deg : float Grid rotation angle in degrees. Defined relative to the co-ordinate system axis and counter clockwise. Returns ------- easting : 2d-array Array containing the values of the easting coordinates of the grid. northing : 2d-array Array containing the values of the northing coordinates of the grid. Notes ----- The ``x`` and ``y`` coordinates are the abscissa and the ordinate of the unrotated grid before the translation, respectively. """ # Define the grid coordinates before the rotation x = np.linspace(0, spacing[1] * (shape[1] - 1), shape[1]) y = np.linspace(0, spacing[0] * (shape[0] - 1), shape[0]) # Compute a meshgrid x, y = np.meshgrid(x, y) # Rotate and shift to get easting and northing rotation_rad = np.radians(rotation_deg) cos, sin = np.cos(rotation_rad), np.sin(rotation_rad) easting = west + x * cos - y * sin northing = south + x * sin + y * cos return easting, northing