grd2geotiff.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Loop3D/geosoft_grid/blob/main/grd2geotiff.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "37bb9833",
      "metadata": {
        "id": "37bb9833"
      },
      "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": {
        "id": "4124ca77"
      },
      "source": [
        "### Code to convert extracted xarray to geotiff format file"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d786c4c5",
      "metadata": {
        "id": "d786c4c5"
      },
      "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": {
        "id": "360df4ac"
      },
      "source": [
        "### Example code to convert uncompressed grid"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f2934ef2",
      "metadata": {
        "id": "f2934ef2"
      },
      "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": {
        "id": "24f4a5a9"
      },
      "source": [
        "### Example code to convert compressed grid"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e53f5bfa",
      "metadata": {
        "id": "e53f5bfa"
      },
      "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": {
        "id": "f8b5f1c3"
      },
      "source": [
        "### Example code to parse XML to get epsg"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "70b37221",
      "metadata": {
        "id": "70b37221"
      },
      "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.9.15"
    },
    "colab": {
      "provenance": [],
      "include_colab_link": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}