Skip to content
Extraits de code Groupes Projets
reclassify_by_interval.ipynb 5,49 ko
Newer Older
  • Learn to ignore specific revisions
  • Boris Nörgaard's avatar
    Boris Nörgaard a validé
    {
     "cells": [
      {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
        "# Reclassify raster by defined intervals"
       ]
      },
      {
       "cell_type": "code",
       "execution_count": 63,
       "metadata": {},
       "outputs": [
        {
         "name": "stdout",
         "output_type": "stream",
         "text": [
          "All libraries successfully imported!\n"
         ]
        }
       ],
       "source": [
        "import rasterio\n",
        "import numpy as np\n",
        "import os\n",
        "\n",
        "print('All libraries successfully imported!')"
       ]
      },
      {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
        "## Set directory"
       ]
      },
      {
       "cell_type": "code",
       "execution_count": 64,
       "metadata": {},
       "outputs": [],
       "source": [
        "computer_path = '/export/miro/ndeffense/LBRAT2104/'\n",
        "grp_letter    = 'X'\n",
        "\n",
        "# Directory for all work files\n",
        "work_path = f'{computer_path}GROUP_{grp_letter}/WORK/'\n",
        "\n",
        "input_file = f'{work_path}NDVI/T31UFS_20200417T104021_NDVI.tif'\n",
        "\n",
        "output_file = f'{input_file[:-4]}_reclassified.tif'"
       ]
      },
      {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
        "## Set parameters"
       ]
      },
      {
       "cell_type": "code",
       "execution_count": 65,
       "metadata": {},
       "outputs": [
        {
         "name": "stdout",
         "output_type": "stream",
         "text": [
          "[-1, -0.5, 0, 0.5, 1]\n"
         ]
        }
       ],
       "source": [
        "nodata_val = -10000\n",
        "\n",
        "# User must defined intervals\n",
        "interval = [-1,-0.5,0,0.5,1]\n",
        "\n",
        "dtype_out = 'int16'\n",
        "\n",
        "print(interval)\n"
       ]
      },
      {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
        "## Reclassify raster by defined intervals\n",
        "\n",
        "Return the indices of the bins to which each value in input array belongs.\n",
        "\n",
        "\n",
        "| `right` |   **order of bins** |  **returned index `i` satisfies** |\n",
        "| --- | --- |  --- |\n",
        "|``False``|increasing | ``bins[i-1] <= x < bins[i]`` |\n",
        "|``True``| increasing | ``bins[i-1] < x <= bins[i]`` |\n",
        "|``False``|  decreasing | ``bins[i-1] > x >= bins[i]`` |\n",
        "|``True``|   decreasing | ``bins[i-1] >= x > bins[i]`` |\n",
        "\n",
        "\n",
        "By default, `right` = False\n",
        "\n",
        "If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is returned as appropriate.\n",
        "\n",
        "https://numpy.org/doc/stable/reference/generated/numpy.digitize.html"
       ]
      },
      {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
        "### Example"
       ]
      },
      {
       "cell_type": "code",
       "execution_count": 66,
       "metadata": {},
       "outputs": [
        {
         "name": "stdout",
         "output_type": "stream",
         "text": [
          "Original continuous values : [ 1.2 10.  12.4 15.5 20. ]\n",
          "Interval : [ 0  5 10 15 20]\n",
          "Reclassified discrete values : \n",
          "[1 2 3 4 4]\n",
          "[1 3 3 4 5]\n"
         ]
        }
       ],
       "source": [
        "x = np.array([1.2, 10.0, 12.4, 15.5, 20.])\n",
        "bins = np.array([0, 5, 10, 15, 20])\n",
        "\n",
        "inds_right_true = np.digitize(x,bins,right=True)\n",
        "inds_right_false = np.digitize(x,bins,right=False)\n",
        "\n",
        "print(f'Original continuous values : {x}')\n",
        "print(f'Interval : {bins}')\n",
        "print('Reclassified discrete values : ')\n",
        "print(inds_right_true)\n",
        "print(inds_right_false)"
       ]
      },
      {
       "cell_type": "code",
       "execution_count": 67,
       "metadata": {},
       "outputs": [
        {
         "name": "stdout",
         "output_type": "stream",
         "text": [
          "data type = int16\n",
          "[[[4 4 4 ... 4 4 4]\n",
          "  [4 4 4 ... 4 4 3]\n",
          "  [4 4 4 ... 4 3 3]\n",
          "  ...\n",
          "  [4 4 4 ... 3 3 3]\n",
          "  [4 4 4 ... 3 3 3]\n",
          "  [4 4 4 ... 3 3 3]]]\n"
         ]
        }
       ],
       "source": [
        "src = rasterio.open(input_file, 'r')\n",
        "im_arr = src.read()\n",
        "\n",
        "# Update the dtype in raster metadata (= profile)\n",
        "profile = src.profile\n",
        "profile.update(dtype = dtype_out)\n",
        "\n",
        "\n",
        "# Replace -10000 by np.nan\n",
        "im_arr[im_arr==nodata_val] = np.nan\n",
        "\n",
        "# Create a mask with all no data value\n",
        "mask = np.isnan(im_arr)\n",
        "\n",
        "# Convert interval into array\n",
        "bins = np.array(interval)\n",
        "\n",
        "# Return the indices of the bins to which each value in input array belongs\n",
        "im_arr_reclass = np.digitize(im_arr, bins, right=False)\n",
        "\n",
        "# Apply mask on reclassified raster\n",
        "im_arr_reclass = np.where(mask, nodata_val, im_arr_reclass)\n",
        "\n",
        "# Change dtype of raster to match the dtype of profile\n",
        "im_arr_reclass = im_arr_reclass.astype(dtype_out)\n",
        "\n",
        "print(f'data type = {im_arr_reclass.dtype}')\n",
        "\n",
        "print(im_arr_reclass)\n",
        "\n",
        "# Write output file\n",
        "dst = rasterio.open(output_file, \"w\", **profile)\n",
        "dst.write(im_arr_reclass)\n",
        "\n",
        "src.close()\n",
        "dst.close()"
       ]
      }
     ],
     "metadata": {
      "interpreter": {
       "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
      },
      "kernelspec": {
       "display_name": "Python 3.6.12 64-bit",
       "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.6.12"
      },
      "orig_nbformat": 4
     },
     "nbformat": 4,
     "nbformat_minor": 2
    }