Index_mosaic_workflow.ipynb

NDVI-mosaic workflow.ipynb

{
  "cells": [
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:17:37.647849Z",
          "end_time": "2022-02-25T07:17:38.171313Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "import os\nimport owslib.wcs as wcs\nimport numpy as np\nimport rasterio as rio\nfrom pathlib import Path",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# 1000x1000px mosaics as example data"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:17:38.173603Z",
          "end_time": "2022-02-25T07:17:38.176787Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "ndindex = 'ndvi'\nwcs_url = 'https://data.nsdc.fmi.fi/geoserver/wcs'\nlayer = f'PTA:s2m_{ndindex}'",
      "execution_count": 2,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:17:38.178767Z",
          "end_time": "2022-02-25T07:17:38.183073Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "ix_path = Path(ndindex)\nos.makedirs(ix_path, exist_ok=True)\ndatapath = ix_path/'testdata'\nyears = [2016,2017,2018,2019,2020,2021]\nos.makedirs(datapath, exist_ok=True)\nfor y in years: os.makedirs(datapath/str(y), exist_ok=True)",
      "execution_count": 3,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:17:38.184826Z",
          "end_time": "2022-02-25T07:17:38.598674Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "w = wcs.WebCoverageService(wcs_url, version='1.0.0')\nc = w[layer]\nt = sorted(c.timepositions)\nt = [ts for ts in t if '15' not in ts]",
      "execution_count": 4,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:17:38.600246Z",
          "end_time": "2022-02-25T07:17:38.602695Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "bounds = (580000,7660000,590000,7670000)",
      "execution_count": 5,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Subset data. If data is already downloaded and monthly no need to download it from API. Here we get only one mosaic per month (first-last), there are also mosaics for middle-of-month intervals if needed."
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:17:38.604580Z",
          "end_time": "2022-02-25T07:18:12.702080Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "size = (int((bounds[2] - bounds[0])/10), int((bounds[3] - bounds[1])/10))\nfor ts in t:\n    outfn = f'{ts.split(\"T\")[0][:-2].replace(\"-\",\"\")}.tif'\n    img = w.getCoverage(layer,\n                   crs='EPSG:3067',\n                   bbox=bounds, \n                   WIDTH=size[0],\n                   HEIGHT=size[1],\n                   format='GeoTIFF',\n                   time=[ts])\n    out = open(datapath/ts[:4]/outfn, 'wb')\n    out.write(img.read())\n    out.close()",
      "execution_count": 6,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Fill gaps with the maximum ndvi of two previous years. This way we can produce maps for 2018, 2019, 2020 and 2021"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:12.703114Z",
          "end_time": "2022-02-25T07:18:12.707497Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "proc_paths = ix_path/'filled'\nos.makedirs(proc_paths, exist_ok=True)\nfillyears = [2018,2019,2020,2021]\nfor f in fillyears: os.makedirs(proc_paths/str(f), exist_ok=True)",
      "execution_count": 7,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:12.710731Z",
          "end_time": "2022-02-25T07:18:13.287230Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "for f in fillyears:\n    for mos in os.listdir(datapath/str(f)):\n        with rio.open(datapath/str(f)/mos) as src:\n            vals = src.read()[0]\n            meta = src.meta\n        with rio.open(datapath/str(f-1)/mos.replace(str(f), str(f-1))) as src_1:\n            vals_1 = src_1.read()[0]\n        with rio.open(datapath/str(f-2)/mos.replace(str(f), str(f-2))) as src_2:\n            vals_2 = src_2.read()[0]\n        vals[vals==0] = np.maximum(vals_1, vals_2)[vals==0]\n        with rio.open(proc_paths/str(f)/mos, 'w', **meta) as dest:\n            dest.write(vals, 1)",
      "execution_count": 8,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Interpolate with `rasterio.fill.fillnodata`, 3 iterations, 100x100px window"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:13.288754Z",
          "end_time": "2022-02-25T07:18:13.292828Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "interp_paths = ix_path/'interp'\nos.makedirs(interp_paths, exist_ok=True)\nfor f in fillyears: os.makedirs(interp_paths/str(f), exist_ok=True)",
      "execution_count": 9,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:13.294680Z",
          "end_time": "2022-02-25T07:18:22.817357Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "from rasterio.fill import fillnodata\n\nfor f in fillyears:\n    for mos in os.listdir(proc_paths/str(f)):\n        with rio.open(proc_paths/str(f)/mos) as src:\n            prof = src.profile\n            arr = src.read(1)\n            arr_filled = fillnodata(arr, mask=src.read_masks(1), max_search_distance=200, smoothing_iterations=3)\n        with rio.open(interp_paths/str(f)/mos, 'w', **prof) as dest:\n            dest.write_band(1, arr_filled)",
      "execution_count": 10,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Fill April, May and October with NDVI from winter mosaic"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:22.819554Z",
          "end_time": "2022-02-25T07:18:23.057339Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "for f in fillyears:\n    with rio.open(f'Winter_median_{ndindex}.tif') as wm:\n        wm_vals = wm.read()[0]\n        for mos in [m for m in os.listdir(interp_paths/str(f)) if any(mon in m for  mon in ('04', '05', '10'))]:\n            with rio.open(interp_paths/str(f)/mos) as src:\n                data = src.read()[0]\n                prof = src.profile\n            data[data==0] = wm_vals[data==0]\n            with rio.open(interp_paths/str(f)/mos, 'w', **prof) as dest:\n                dest.write_band(1, data)\n            ",
      "execution_count": 11,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Fill remaining gaps in August with the average of July and September of the same year"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:23.059103Z",
          "end_time": "2022-02-25T07:18:23.210366Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "for f in fillyears:\n    with rio.open(interp_paths/str(f)/f'{f}08.tif') as src:\n        prof = src.profile\n        aug = src.read()[0]\n    with rio.open(interp_paths/str(f)/f'{f}07.tif') as src:\n        jul = src.read()[0]\n    with rio.open(interp_paths/str(f)/f'{f}09.tif') as src:\n        sep = src.read()[0]\n    aug[aug==0] = np.mean(np.array([jul, sep]), axis=0)[aug==0]\n    with rio.open(interp_paths/str(f)/f'{f}08.tif', 'w', **prof) as dest:\n        dest.write_band(1, aug)",
      "execution_count": 12,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Make statistics rasters (yearly min, max, mean and median)."
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:23.212223Z",
          "end_time": "2022-02-25T07:18:23.217196Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "statspath = ix_path/'stats'\nos.makedirs(statspath, exist_ok=True)\nfor f in fillyears: os.makedirs(statspath/str(f), exist_ok=True)",
      "execution_count": 13,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Compute yearly stats:\n* `mean`\n* `median`\n* `max` \n* `min`\n* `stdev`\n* `sum` as a proxy for total productivity"
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:23.219973Z",
          "end_time": "2022-02-25T07:18:24.542064Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "for f in fillyears:\n    mosaics = []\n    for mos in os.listdir(interp_paths/str(f)):\n        with rio.open(interp_paths/str(f)/mos) as src:\n            prof = src.profile\n            data = src.read()[0]\n        mosaics.append(data)\n    mosaics = np.array(mosaics)\n    with rio.open(statspath/str(f)/'mean.tif', 'w', **prof) as dest:\n        dest.write(mosaics.mean(axis=0),1)\n    with rio.open(statspath/str(f)/'median.tif', 'w', **prof) as dest:\n        dest.write(np.median(mosaics, axis=0),1)\n    with rio.open(statspath/str(f)/'min.tif', 'w', **prof) as dest:\n        dest.write(mosaics.min(axis=0),1)\n    with rio.open(statspath/str(f)/'max.tif', 'w', **prof) as dest:\n        dest.write(mosaics.max(axis=0),1)\n    prof.update({'dtype': 'float32',\n                 'nodataval': None})\n    with rio.open(statspath/str(f)/'stdev.tif', 'w', **prof) as dest:\n        dest.write(mosaics.astype(np.float32).std(axis=0), 1)\n    prof.update({'dtype':'uint16',\n                 'nodataval': 0})\n    mosaics = mosaics.astype(np.int16)\n    with rio.open(statspath/str(f)/'sum.tif', 'w', **prof) as dest:\n        dest.write(mosaics.sum(axis=0), 1)\n    prof.update({'dtype':'int16',\n                 'nodataval': 0})\n    with rio.open(statspath/str(f)/'amp.tif', 'w', **prof) as dest:\n        dest.write(mosaics.max(axis=0)-np.median(mosaics, axis=0), 1)",
      "execution_count": 14,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Monthly amplitudes NDVI$_{max}$ - NDVI$_{base}$, where NDVI$_{base}$ is the pixelwise median of all input data layers. We need to set nodatavalue to `None` and datatype to `int16` as the values can go below zero."
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:24.543800Z",
          "end_time": "2022-02-25T07:18:24.548403Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "amp_path = ix_path/'amp'\nos.makedirs(amp_path, exist_ok=True)\nfor f in fillyears: os.makedirs(amp_path/str(f), exist_ok=True)",
      "execution_count": 15,
      "outputs": []
    },
    {
      "metadata": {
        "ExecuteTime": {
          "start_time": "2022-02-25T07:18:24.550309Z",
          "end_time": "2022-02-25T07:18:24.907483Z"
        },
        "trusted": true
      },
      "cell_type": "code",
      "source": "for f in fillyears:\n    with rio.open(statspath/str(f)/'median.tif') as med:\n        medvals = med.read()[0]\n        medvals = medvals.astype(np.int16)\n        for mos in os.listdir(interp_paths/str(f)):\n            with rio.open(interp_paths/str(f)/mos) as src:\n                vals = src.read()[0]\n                prof = src.profile\n            prof.update({'dtype': 'int16',\n                         'nodata': None})\n            vals = vals.astype(np.int16)\n            amp = vals - medvals\n            with rio.open(amp_path/str(f)/mos, 'w', **prof) as dest:\n                dest.write(amp, 1)",
      "execution_count": 16,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "toc": {
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": false,
      "base_numbering": 1,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {},
      "toc_section_display": true,
      "toc_window_display": false
    },
    "latex_envs": {
      "eqNumInitial": 1,
      "eqLabelWithNumbers": true,
      "current_citInitial": 1,
      "cite_by": "apalike",
      "bibliofile": "biblio.bib",
      "LaTeX_envs_menu_present": true,
      "labels_anchors": false,
      "latex_user_defs": false,
      "user_envs_cfg": false,
      "report_style_numbering": false,
      "autoclose": false,
      "autocomplete": true,
      "hotkeys": {
        "equation": "Ctrl-E",
        "itemize": "Ctrl-I"
      }
    },
    "language_info": {
      "name": "python",
      "version": "3.9.9",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "gist": {
      "id": "0a27163fdd4d1e7b8850fa2308bfaa2e",
      "data": {
        "description": "Index_mosaic_workflow.ipynb",
        "public": true
      }
    },
    "_draft": {
      "nbviewer_url": "https://gist.github.com/0a27163fdd4d1e7b8850fa2308bfaa2e"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}