From Fatianda (use_this).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37bb9833",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright (c) 2018 The Harmonica Developers.\n",
    "# Distributed under the terms of the BSD 3-Clause License.\n",
    "# SPDX-License-Identifier: BSD-3-Clause\n",
    "#\n",
    "# This code is part of the Fatiando a Terra project (https://www.fatiando.org)\n",
    "#\n",
    "# The support for compressed grids was inspired by the code in\n",
    "# https://github.com/Loop3D/geosoft_grid copyrighted by Loop3D and released\n",
    "# under the MIT License.\n",
    "\"\"\"\n",
    "Function to read Oasis Montaj© .grd file\n",
    "\"\"\"\n",
    "\n",
    "import array\n",
    "import zlib\n",
    "\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "# Define the valid element sizes (ES variable) for GRD files\n",
    "# (values > 1024 correspond to compressed versions of the grid)\n",
    "VALID_ELEMENT_SIZES = (1, 2, 4, 8, 1024 + 1, 1024 + 2, 1024 + 4, 1024 + 8)\n",
    "# Define dummy values for each data type\n",
    "DUMMIES = {\n",
    "    \"b\": -127,\n",
    "    \"B\": 255,\n",
    "    \"h\": -32767,\n",
    "    \"H\": 65535,\n",
    "    \"i\": -2147483647,\n",
    "    \"I\": 4294967295,\n",
    "    \"f\": -1e32,\n",
    "    \"d\": -1e32,\n",
    "}\n",
    "\n",
    "\n",
    "def load_oasis_montaj_grid(fname):\n",
    "    \"\"\"\n",
    "    Reads gridded data from an Oasis Montaj© .grd file.\n",
    "    The version 2 of the Geosoft© Grid File Format (GRD) stores gridded\n",
    "    products in binary data. This function can read those files and parse the\n",
    "    information in the header. It returns the data in\n",
    "    a :class:`xarray.DataArray` for convenience.\n",
    "    .. warning::\n",
    "        This function has not been tested against a wide range of GRD files.\n",
    "        This could lead to incorrect readings of the stored data. Please report\n",
    "        any unwanted behaviour by opening an issue in Harmonica:\n",
    "        https://github.com/fatiando/harmonica/issues\n",
    "    .. important::\n",
    "        This function only supports reading GRD files using the **version 2**\n",
    "        of the Geosoft© Grid File Format.\n",
    "    .. important::\n",
    "        This function is not supporting orderings different than ±1,\n",
    "        or colour grids.\n",
    "    Parameters\n",
    "    ----------\n",
    "    fname : string or file-like object\n",
    "        Path to the .grd file.\n",
    "    Returns\n",
    "    -------\n",
    "    grid : :class:`xarray.DataArray`\n",
    "        :class:`xarray.DataArray` containing the grid, its coordinates and\n",
    "        header information.\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    https://help.seequent.com/Oasis-montaj/9.9/en/Content/ss/glossary/grid_file_format__grd.htm\n",
    "    https://github.com/Loop3D/geosoft_grid\n",
    "    \"\"\"\n",
    "    # Read the header and the grid array\n",
    "    with open(fname, \"rb\") as grd_file:\n",
    "        # Read the header (first 512 bytes)\n",
    "        header = _read_header(grd_file.read(512))\n",
    "        # Check for valid flags\n",
    "        _check_ordering(header[\"ordering\"])\n",
    "        _check_sign_flag(header[\"sign_flag\"])\n",
    "        # Get data type for the grid elements\n",
    "        data_type = _get_data_type(header[\"n_bytes_per_element\"], header[\"sign_flag\"])\n",
    "        # Read grid\n",
    "        grid = grd_file.read()\n",
    "    # Decompress grid if needed\n",
    "    if header[\"n_bytes_per_element\"] > 1024:\n",
    "        grid = _decompress_grid(grid)\n",
    "    # Load the grid values as an array with the proper data_type\n",
    "    grid = array.array(data_type, grid)\n",
    "    # Convert to numpy array as float64\n",
    "    grid = np.array(grid, dtype=np.float64)\n",
    "    # Remove dummy values\n",
    "    grid = _remove_dummies(grid, data_type)\n",
    "    # Scale the grid\n",
    "    grid = np.array(grid / header[\"data_factor\"] + header[\"base_value\"])\n",
    "    # Reshape the grid based on the ordering\n",
    "    if header[\"ordering\"] == 1:\n",
    "        order = \"C\"\n",
    "        shape = (header[\"shape_v\"], header[\"shape_e\"])\n",
    "        spacing = (header[\"spacing_v\"], header[\"spacing_e\"])\n",
    "    elif header[\"ordering\"] == -1:\n",
    "        order = \"F\"\n",
    "        shape = (header[\"shape_e\"], header[\"shape_v\"])\n",
    "        spacing = (header[\"spacing_e\"], header[\"spacing_v\"])\n",
    "    grid = grid.reshape(shape, order=order)\n",
    "    # Build coords\n",
    "    if header[\"rotation\"] == 0:\n",
    "        easting, northing = _build_coordinates(\n",
    "            header[\"x_origin\"], header[\"y_origin\"], shape, spacing\n",
    "        )\n",
    "        dims = (\"northing\", \"easting\")\n",
    "        coords = {\"easting\": easting, \"northing\": northing}\n",
    "    else:\n",
    "        easting, northing = _build_rotated_coordinates(\n",
    "            header[\"x_origin\"], header[\"y_origin\"], shape, spacing, header[\"rotation\"]\n",
    "        )\n",
    "        dims = (\"y\", \"x\")\n",
    "        coords = {\"easting\": (dims, easting), \"northing\": (dims, northing)}\n",
    "    # Build an xarray.DataArray for the grid\n",
    "    grid = xr.DataArray(\n",
    "        grid,\n",
    "        coords=coords,\n",
    "        dims=dims,\n",
    "        attrs=header,\n",
    "    )\n",
    "    return grid\n",
    "\n",
    "\n",
    "def _read_header(header_bytes):\n",
    "    \"\"\"\n",
    "    Read GRD file header\n",
    "    Parameters\n",
    "    ----------\n",
    "    header_bytes : byte\n",
    "        A sequence of 512 bytes containing the header of a\n",
    "        GRD file.\n",
    "    Returns\n",
    "    -------\n",
    "    header : dict\n",
    "        Dictionary containing the information present in the\n",
    "        header.\n",
    "    Notes\n",
    "    -----\n",
    "    The GRD header consists in 512 contiguous bytes.\n",
    "    It's divided in four sections:\n",
    "    * Data Storage\n",
    "    * Geographic Information\n",
    "    * Data (Z) Scaling\n",
    "    * Undefined Application Parameters\n",
    "    \"\"\"\n",
    "    header = {}\n",
    "    # Read data storage\n",
    "    ES, SF, NE, NV, KX = array.array(\"i\", header_bytes[0 : 5 * 4])  # noqa: N806\n",
    "    header.update(\n",
    "        {\n",
    "            \"n_bytes_per_element\": ES,\n",
    "            \"sign_flag\": SF,\n",
    "            \"shape_e\": NE,\n",
    "            \"shape_v\": NV,\n",
    "            \"ordering\": KX,\n",
    "        }\n",
    "    )\n",
    "    # Read geographic info\n",
    "    DE, DV, X0, Y0, ROT = array.array(\"d\", header_bytes[20 : 20 + 5 * 8])  # noqa: N806\n",
    "    header.update(\n",
    "        {\n",
    "            \"spacing_e\": DE,\n",
    "            \"spacing_v\": DV,\n",
    "            \"x_origin\": X0,\n",
    "            \"y_origin\": Y0,\n",
    "            \"rotation\": ROT,\n",
    "        }\n",
    "    )\n",
    "    # Read data scaling\n",
    "    ZBASE, ZMULT = array.array(\"d\", header_bytes[60 : 60 + 2 * 8])  # noqa: N806\n",
    "    header.update(\n",
    "        {\n",
    "            \"base_value\": ZBASE,\n",
    "            \"data_factor\": ZMULT,\n",
    "        }\n",
    "    )\n",
    "    # Read optional parameters\n",
    "    # (ignore map LABEL and MAPNO)\n",
    "    PROJ, UNITX, UNITY, UNITZ, NVPTS = array.array(  # noqa: N806\n",
    "        \"i\", header_bytes[140 : 140 + 5 * 4]\n",
    "    )\n",
    "    IZMIN, IZMAX, IZMED, IZMEA = array.array(  # noqa: N806\n",
    "        \"f\", header_bytes[160 : 160 + 4 * 4]\n",
    "    )\n",
    "    (ZVAR,) = array.array(\"d\", header_bytes[176 : 176 + 8])  # noqa: N806\n",
    "    (PRCS,) = array.array(\"i\", header_bytes[184 : 184 + 4])  # noqa: N806\n",
    "    header.update(\n",
    "        {\n",
    "            \"map_projection\": PROJ,\n",
    "            \"units_x\": UNITX,\n",
    "            \"units_y\": UNITY,\n",
    "            \"units_z\": UNITZ,\n",
    "            \"n_valid_points\": NVPTS,\n",
    "            \"grid_min\": IZMIN,\n",
    "            \"grid_max\": IZMAX,\n",
    "            \"grid_median\": IZMED,\n",
    "            \"grid_mean\": IZMEA,\n",
    "            \"grid_variance\": ZVAR,\n",
    "            \"process_flag\": PRCS,\n",
    "        }\n",
    "    )\n",
    "    return header\n",
    "\n",
    "\n",
    "def _check_ordering(ordering):\n",
    "    \"\"\"\n",
    "    Check if the ordering value is within the ones we are supporting\n",
    "    \"\"\"\n",
    "    if ordering not in (-1, 1):\n",
    "        raise NotImplementedError(\n",
    "            f\"Found an ordering (a.k.a as KX) equal to '{ordering}'. \"\n",
    "            + \"Only orderings equal to 1 and -1 are supported.\"\n",
    "        )\n",
    "\n",
    "\n",
    "def _check_sign_flag(sign_flag):\n",
    "    \"\"\"\n",
    "    Check if sign_flag value is within the ones we are supporting\n",
    "    \"\"\"\n",
    "    if sign_flag == 3:\n",
    "        raise NotImplementedError(\n",
    "            \"Reading .grd files with colour grids is not currenty supported.\"\n",
    "        )\n",
    "\n",
    "\n",
    "def _get_data_type(n_bytes_per_element, sign_flag):\n",
    "    \"\"\"\n",
    "    Return the data type for the grid values\n",
    "    References\n",
    "    ----------\n",
    "    https://docs.python.org/3/library/array.html\n",
    "    \"\"\"\n",
    "    # Check if number of bytes per element is valid\n",
    "    if n_bytes_per_element not in VALID_ELEMENT_SIZES:\n",
    "        raise NotImplementedError(\n",
    "            \"Found a 'Grid data element size' (a.k.a. 'ES') value \"\n",
    "            + f\"of '{n_bytes_per_element}'. \"\n",
    "            \"Only values equal to 1, 2, 4 and 8 are valid, \"\n",
    "            + \"along with their compressed counterparts (1025, 1026, 1028, 1032).\"\n",
    "        )\n",
    "    # Shift the n_bytes_per_element in case of compressed grids\n",
    "    if n_bytes_per_element > 1024:\n",
    "        n_bytes_per_element -= 1024\n",
    "    # Determine the data type of the grid elements\n",
    "    if n_bytes_per_element == 1:\n",
    "        if sign_flag == 0:\n",
    "            data_type = \"B\"  # unsigned char\n",
    "        elif sign_flag == 1:\n",
    "            data_type = \"b\"  # signed char\n",
    "    elif n_bytes_per_element == 2:\n",
    "        if sign_flag == 0:\n",
    "            data_type = \"H\"  # unsigned short\n",
    "        elif sign_flag == 1:\n",
    "            data_type = \"h\"  # signed short\n",
    "    elif n_bytes_per_element == 4:\n",
    "        if sign_flag == 0:\n",
    "            data_type = \"I\"  # unsigned int\n",
    "        elif sign_flag == 1:\n",
    "            data_type = \"i\"  # signed int\n",
    "        elif sign_flag == 2:\n",
    "            data_type = \"f\"  # float\n",
    "    elif n_bytes_per_element == 8:\n",
    "        data_type = \"d\"\n",
    "    return data_type\n",
    "\n",
    "\n",
    "def _remove_dummies(grid, data_type):\n",
    "    \"\"\"\n",
    "    Replace dummy values for NaNs\n",
    "    \"\"\"\n",
    "    # Create dictionary with dummy value for each data type\n",
    "    if data_type in (\"f\", \"d\"):\n",
    "        grid[grid <= DUMMIES[data_type]] = np.nan\n",
    "        return grid\n",
    "    grid[grid == DUMMIES[data_type]] = np.nan\n",
    "    return grid\n",
    "\n",
    "\n",
    "def _decompress_grid(grid_compressed):\n",
    "    \"\"\"\n",
    "    Decompress the grid using gzip\n",
    "    Even if the header specifies that the grid is compressed using a LZRW1\n",
    "    algorithm, it's using gzip instead. The first two 4 bytes sequences\n",
    "    correspond to the compression signature and\n",
    "    to the compression type. We are going to ignore those and start reading\n",
    "    from the number of blocks (offset 8).\n",
    "    Parameters\n",
    "    ----------\n",
    "    grid_compressed : bytes\n",
    "        Sequence of bytes corresponding to the compressed grid. They should be\n",
    "        every byte starting from offset 512 of the GRD file until its end.\n",
    "    Returns\n",
    "    -------\n",
    "    grid : bytes\n",
    "        Uncompressed version of the ``grid_compressed`` parameter.\n",
    "    \"\"\"\n",
    "    # Number of blocks\n",
    "    (n_blocks,) = array.array(\"i\", grid_compressed[8 : 8 + 4])\n",
    "    # Number of vectors per block\n",
    "    (vectors_per_block,) = array.array(\"i\", grid_compressed[12 : 12 + 4])\n",
    "    # File offset from start of every block\n",
    "    block_offsets = array.array(\"q\", grid_compressed[16 : 16 + n_blocks * 8])\n",
    "    # Compressed size of every block\n",
    "    compressed_block_sizes = array.array(\n",
    "        \"i\",\n",
    "        grid_compressed[16 + n_blocks * 8 : 16 + n_blocks * 8 + n_blocks * 4],\n",
    "    )\n",
    "    # Combine grid\n",
    "    grid = b\"\"\n",
    "    # Read each block\n",
    "    for i in range(n_blocks):\n",
    "        # Define the start and end offsets for each compressed blocks\n",
    "        # We need to remove the 512 to account for the missing header.\n",
    "        # There is an unexplained 16 byte header that we also need to remove.\n",
    "        start_offset = block_offsets[i] - 512 + 16\n",
    "        end_offset = compressed_block_sizes[i] + block_offsets[i] - 512\n",
    "        # Decompress the block\n",
    "        grid_sub = zlib.decompress(\n",
    "            grid_compressed[start_offset:end_offset],\n",
    "            bufsize=zlib.DEF_BUF_SIZE,\n",
    "        )\n",
    "        # Add it to the running grid\n",
    "        grid += grid_sub\n",
    "    return grid\n",
    "\n",
    "\n",
    "def _build_coordinates(west, south, shape, spacing):\n",
    "    \"\"\"\n",
    "    Create the coordinates for the grid\n",
    "    Generates 1d arrays for the easting and northing coordinates of the grid.\n",
    "    Assumes unrotated grids.\n",
    "    Parameters\n",
    "    ----------\n",
    "    west : float\n",
    "        Westernmost coordinate of the grid.\n",
    "    south : float\n",
    "        Southernmost coordinate of the grid.\n",
    "    shape : tuple\n",
    "        Tuple of ints containing the number of elements along each direction in\n",
    "        the following order: ``n_northing``, ``n_easting``\n",
    "    spacing : tuple\n",
    "        Tuple of floats containing the distance between adjacent grid elements\n",
    "        along each direction in the following order: ``s_northing``,\n",
    "        ``s_easting``.\n",
    "    Returns\n",
    "    -------\n",
    "    easting : 1d-array\n",
    "        Array containing the values of the easting coordinates of the grid.\n",
    "    northing : 1d-array\n",
    "        Array containing the values of the northing coordinates of the grid.\n",
    "    \"\"\"\n",
    "    easting = np.linspace(west, west + spacing[1] * (shape[1] - 1), shape[1])\n",
    "    northing = np.linspace(south, south + spacing[0] * (shape[0] - 1), shape[0])\n",
    "    return easting, northing\n",
    "\n",
    "\n",
    "def _build_rotated_coordinates(west, south, shape, spacing, rotation_deg):\n",
    "    \"\"\"\n",
    "    Create the coordinates for a rotated grid\n",
    "    Generates 2d arrays for the easting and northing coordinates of the grid.\n",
    "    Assumes rotated grids.\n",
    "    Parameters\n",
    "    ----------\n",
    "    west : float\n",
    "        Westernmost coordinate of the grid.\n",
    "    south : float\n",
    "        Southernmost coordinate of the grid.\n",
    "    shape : tuple\n",
    "        Tuple of ints containing the number of elements along each unrotated\n",
    "        direction in the following order: ``n_y``, ``n_x``\n",
    "    spacing : tuple\n",
    "        Tuple of floats containing the distance between adjacent grid elements\n",
    "        along each unrotated direction in the following order: ``spacing_y``,\n",
    "        ``spacing_x``.\n",
    "    Returns\n",
    "    -------\n",
    "    easting : 2d-array\n",
    "        Array containing the values of the easting coordinates of the grid.\n",
    "    northing : 2d-array\n",
    "        Array containing the values of the northing coordinates of the grid.\n",
    "    Notes\n",
    "    -----\n",
    "    The ``x`` and ``y`` coordinates are the abscissa and the ordinate of the\n",
    "    unrotated grid before the translation, respectively.\n",
    "    \"\"\"\n",
    "    # Define the grid coordinates before the rotation\n",
    "    x = np.linspace(0, spacing[1] * (shape[1] - 1), shape[1])\n",
    "    y = np.linspace(0, spacing[0] * (shape[0] - 1), shape[0])\n",
    "    # Compute a meshgrid\n",
    "    x, y = np.meshgrid(x, y)\n",
    "    # Rotate and shift to get easting and northing\n",
    "    rotation_rad = np.radians(rotation_deg)\n",
    "    cos, sin = np.cos(rotation_rad), np.sin(rotation_rad)\n",
    "    easting = west + x * cos - y * sin\n",
    "    northing = south + x * sin + y * cos\n",
    "    return easting, northing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4124ca77",
   "metadata": {},
   "source": [
    "### Code to convert extracted xarray to geotiff format file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d786c4c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterio.transform import from_origin\n",
    "\n",
    "def write_geotiff(out_path,grid_uc,epsg):\n",
    "    \"\"\"\n",
    "    write out grid in geotiff format\n",
    "    Parameters\n",
    "    ----------\n",
    "    out_path : string\n",
    "        file path and name to write to\n",
    "    grid_uc : xarray\n",
    "        grid and header info  \n",
    "    epsg: string\n",
    "        user defined epsg\n",
    "    \"\"\"    \n",
    "    transform = from_origin(grid_uc.attrs['x_origin']-(grid_uc.attrs['spacing_e']/2), grid_uc.attrs['y_origin']-(grid_uc.attrs['spacing_v']/2), grid_uc.attrs['spacing_e'], -grid_uc.attrs['spacing_v'])\n",
    "\n",
    "    new_dataset = rasterio.open(\n",
    "        out_path,\n",
    "        \"w\",\n",
    "        driver=\"GTiff\",\n",
    "        height=grid_uc.attrs['shape_v'],\n",
    "        width=grid_uc.attrs['shape_e'],\n",
    "        count=1,\n",
    "        dtype=rasterio.float32,\n",
    "        crs='epsg:'+epsg,#\"+proj=longlat\",\n",
    "        transform=transform,\n",
    "        nodata=-1e8,\n",
    "    )\n",
    "\n",
    "    new_dataset.write(grid_uc.values, 1)\n",
    "    new_dataset.close()\n",
    "    print(\"dtm geotif saved as\", out_path)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "360df4ac",
   "metadata": {},
   "source": [
    "### Example code to convert uncompressed grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2934ef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename='./test_data/test_uncompressed.grd' # uncompressed example\n",
    "grid_uc=load_oasis_montaj_grid(filename)\n",
    "write_geotiff('./test1.tif',grid_uc,'4326')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24f4a5a9",
   "metadata": {},
   "source": [
    "### Example code to convert compressed grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e53f5bfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename='./test_data/test_compressed.grd' # compressed example\n",
    "grid_c=load_oasis_montaj_grid(filename)\n",
    "write_geotiff('./test2.tif',grid_c,'4326')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8b5f1c3",
   "metadata": {},
   "source": [
    "### Example code to parse XML to get epsg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70b37221",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Suggested code to parse xml to get CRS from\n",
    "# Santiago Soler \n",
    "# https://github.com/fatiando/harmonica/pull/348#issuecomment-1327811755\n",
    "\n",
    "\n",
    "def extract_proj_str(fname):\n",
    "    with open(fname, \"r\") as f:\n",
    "        for line in f:\n",
    "            if \"projection=\" in line:\n",
    "                # print lines with projection info\n",
    "                print(line)\n",
    "\n",
    "                # extract projection string\n",
    "                proj = line.split('projection=\"')[1].split('\" ')[0]\n",
    "                print(proj)\n",
    "\n",
    "                # remove non-alphanumeric characters if present\n",
    "                if proj.isalnum() is False:\n",
    "                    proj = ''.join(filter(str.isalnum, proj))\n",
    "\n",
    "                assert proj.isalnum\n",
    "    try:\n",
    "        proj\n",
    "    except:\n",
    "        raise NameError(\"string 'projection=' not found in file.\")\n",
    "        proj = None\n",
    "\n",
    "    return proj\n",
    "\n",
    "proj = extract_proj_str(filename+\".xml\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}