TEAK_geo3_giannico_method_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook quantizes the footprints into equal intervals of 5% based on the footprint cell weights.  Before quantizing, it upscales to the resolution of, and aligns pixels with the hyperspectral data. After Quantizing it records all intervals up to the 0.95 quantile as separate boolean masks in a stacked data array. Layer 0 corresponds to the 0.05 - 0.10 quantile contribution zone, layer 19 corresponds to 0.95 - 1, i.e. the cells in layer 19 are those that are at or above the 95th percentile of contributors to the flux at the eddy tower.   \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import rioxarray\n",
    "from xrspatial.classify import equal_interval\n",
    "import os\n",
    "from tqdm import tqdm \n",
    "import matplotlib.pyplot as plt\n",
    "from rasterio import uint8, features\n",
    "import json\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Polygon, mapping\n",
    "import fiona\n",
    "from fiona.crs import from_epsg\n",
    "import rasterio\n",
    "from pandas import to_datetime\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "\n",
    "# paths TODO: use os.path.join to create\n",
    "site = 'TEAK'\n",
    "basepath = f'/media/data/NEON/{site}'\n",
    "ncdf = '/media/data/NEON/TEAK/mosaic/TEAK_2019_mosaic.nc'\n",
    "foot_path = '/media/data/NEON/TEAK/footprints'\n",
    "\n",
    "test = os.path.join(basepath, 'test.tiff')\n",
    "\n",
    "# get names of footprints\n",
    "footprints = [os.path.join(foot_path, f) for f in os.listdir(foot_path) if 'summary' not in f and '.xml' not in f]\n",
    "\n",
    "# open netcdf mosaic\n",
    "data = xr.open_dataset(ncdf)\n",
    "\n",
    "# set epsg from attrs\n",
    "data = data.rio.write_crs(data.attrs['epsg'])\n",
    "\n",
    "# clip the hyperspectral derived data to footprint extent\n",
    "footprint = rioxarray.open_rasterio(footprints[0])\n",
    "data = data.rio.clip_box(minx=np.atleast_1d(footprint.x.min())[0],\n",
    "                        miny=np.atleast_1d(footprint.y.min())[0],\n",
    "                        maxx=np.atleast_1d(footprint.x.max())[0],\n",
    "                        maxy=np.atleast_1d(footprint.y.max())[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ensure directory exists for vectors to live in\n",
    "footmask_home = os.path.join(basepath, 'quantile_footprints')\n",
    "os.makedirs(footmask_home, exist_ok=True)\n",
    "\n",
    "\n",
    "\n",
    "for foot in tqdm(footprints[:1]):\n",
    "    # get timestamp\n",
    "    t = to_datetime(os.path.basename(foot).split('.')[0])\n",
    "\n",
    "    # open footprint\n",
    "    footprint = rioxarray.open_rasterio(foot)\n",
    "\n",
    "    # footprint.shape is(1, y, x) so bust it out 3rd dim\n",
    "    footprint = footprint[0]\n",
    "\n",
    "    # rotate it because it is read in 90 off for some reason\n",
    "    footprint.data = np.rot90(footprint.data, k=3)\n",
    "\n",
    "    # then upscale\n",
    "    footprint = footprint.interp_like(data.kmeans_label)\n",
    "\n",
    "    #  transpose dims\n",
    "    #footprint = footprint.transpose('x', 'y')\n",
    "\n",
    "    # renormalize , interp messes up the sum\n",
    "    footprint.data = footprint.data / np.nansum(footprint.data)\n",
    "\n",
    "    # create equal interval contour classifcations\n",
    "    footprint = equal_interval(footprint, k=20)\n",
    "\n",
    "    # use the equal interval contours to calculate the \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 26%|██▌       | 955/3672 [1:58:48<4:52:45,  6.47s/it] "
     ]
    }
   ],
   "source": [
    "def interp_foot_2_data(footprint, data):\n",
    "   \n",
    "    # footprints are different resolution so we must upscale, first get dims\n",
    "    old_x = footprint.shape[1]\n",
    "    old_y = footprint.shape[0]\n",
    "\n",
    "    # then upscale\n",
    "    footprint = footprint.interp_like(data.kmeans_label)\n",
    "\n",
    "    # fill nulls\n",
    "    footprint = footprint.fillna(0)\n",
    "\n",
    "    # renormalize , interp messes up the sum\n",
    "    footprint = footprint / np.nansum(footprint)\n",
    "\n",
    "    return(footprint)\n",
    "\n",
    "\n",
    "def square_crop(q):\n",
    "    '''returns 1km square datarray cropped around the extent of non-zero entries'''\n",
    "    # find the extent of the q95 mask\n",
    "    xmin = q.x.values[np.argwhere(q.data > 0)[:,1].min()]\n",
    "    xmax = q.x.values[np.argwhere(q.data > 0)[:,1].max()]\n",
    "\n",
    "    ymin = q.y.values[np.argwhere(q.data > 0)[:,0].min()]\n",
    "    ymax = q.y.values[np.argwhere(q.data > 0)[:,0].max()]\n",
    "\n",
    "    # buffer it out to km^2\n",
    "    xbuf = (1000 - (xmax - xmin)) / 2\n",
    "    xmin = xmin - xbuf\n",
    "    xmax = xmax + xbuf\n",
    "\n",
    "    ybuf = (1000 - (ymax - ymin)) / 2\n",
    "    ymin = ymin - ybuf\n",
    "    ymax = ymax + ybuf\n",
    "\n",
    "\n",
    "    mask_x = (q.x >= xmin) & (q.x <= xmax)\n",
    "    mask_y = (q.y >= ymin) & (q.y <= ymax)\n",
    "\n",
    "    crop = q.where(mask_x & mask_y, drop=True)\n",
    "    return crop\n",
    "\n",
    "\n",
    "def initiate_footstack(foot, data):\n",
    "    '''\n",
    "    uses the first footprint to create a new dataset\n",
    "    '''\n",
    "    # get timestamp\n",
    "    t = os.path.basename(foot).split('.')[0]\n",
    "\n",
    "    # open footprint\n",
    "    footprint = rioxarray.open_rasterio(foot)\n",
    "\n",
    "    # footprint.shape is(1, y, x) so bust it out 3rd dim\n",
    "    footprint = footprint[0]\n",
    "\n",
    "    # upscale and align footprint\n",
    "    footprint = interp_foot_2_data(footprint, data)\n",
    "\n",
    "    # create equal interval contour classifcations\n",
    "    footprint = equal_interval(footprint, k=20)\n",
    "\n",
    "    # get the q95 footprint \n",
    "    vals = set(footprint.data[footprint.data > 0].flatten())\n",
    "    val_max = max(vals)\n",
    "    q95 = np.zeros_like(footprint.data, dtype=np.int8)\n",
    "    q95[footprint.data == val_max] = 1\n",
    "    q95 = xr.DataArray(q95, dims=['y', 'x'], coords=[footprint.y, footprint.x])\n",
    "\n",
    "    # crop q95 footprint to 1km square\n",
    "    q95 = square_crop(q95)\n",
    "\n",
    "    # align footprint to q95\n",
    "    global CROP_TEMPLATE\n",
    "    CROP_TEMPLATE, footprint = xr.align(q95, footprint)\n",
    "\n",
    "    # create the dataset\n",
    "    footstack = footprint.to_dataset()\n",
    "    footstack = footstack.drop('equal_interval')\n",
    "    \n",
    "    # ad a new dimension, t, fot timeseries of footrasters \n",
    "    footstack = footstack.expand_dims(dim='t')\n",
    "\n",
    "    # make all 0.5 t0 0.95 quantile rasters\n",
    "    qs = []\n",
    "    names = []\n",
    "\n",
    "    for val in vals:\n",
    "        q = np.zeros_like(footprint.data, dtype=np.int8)\n",
    "        q[footprint.data == val] =  1\n",
    "        q = xr.DataArray(q, dims=['y', 'x'], coords=[footprint.y, footprint.x])\n",
    "        names.append(f'q{int(val * 5)}')\n",
    "        qs.append(q)\n",
    "\n",
    "    foot_series = xr.concat(qs, dim='quantile')\n",
    "\n",
    "    return footstack, foot_series, t\n",
    "\n",
    "\n",
    "def make_foot_series(foot):\n",
    "    \n",
    "    # get timestamp\n",
    "    t = os.path.basename(foot).split('.')[0]\n",
    "\n",
    "    # open footprint\n",
    "    footprint = rioxarray.open_rasterio(foot)\n",
    "\n",
    "    # footprint.shape is(1, y, x) so bust it out 3rd dim\n",
    "    footprint = footprint[0]\n",
    "\n",
    "    # upscale and align footprint\n",
    "    footprint = interp_foot_2_data(footprint, data)\n",
    "\n",
    "    # create equal interval contour classifcations\n",
    "    footprint = equal_interval(footprint, k=20)\n",
    "\n",
    "    # make all 0.5 t0 0.95 quantile rasters\n",
    "    qs = []\n",
    "    names = []\n",
    "\n",
    "    vals = set(footprint.data[footprint.data > 0].flatten())\n",
    "    for val in vals:\n",
    "        q = np.zeros_like(footprint.data, dtype=np.int8)\n",
    "        q[footprint.data == val] =  1\n",
    "        q = xr.DataArray(q, dims=['y', 'x'], coords=[footprint.y, footprint.x])\n",
    "        _ , q = xr.align(CROP_TEMPLATE, q)\n",
    "        names.append(f'q{int(val * 5)}')\n",
    "        qs.append(q)\n",
    "\n",
    "    foot_series = xr.concat(qs, dim='quantile')\n",
    "\n",
    "    return foot_series, t\n",
    "\n",
    "\n",
    "\n",
    "first = True\n",
    "\n",
    "for foot in tqdm(footprints):\n",
    "    if first:\n",
    "        footstack, foot_series, t = initiate_footstack(foot, data)\n",
    "        first = False\n",
    "    else:\n",
    "        foot_series, t = make_foot_series(foot)\n",
    "    \n",
    "    foot_series = foot_series.rename('footprint') \n",
    "    foot_series.to_netcdf(os.path.join(footmask_home, f'{t}.nc'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "0b20d2e798ab320dcd9e3022b3f23321d95437ce0beeb6a34c6f05479d8e44a2"
  },
  "kernelspec": {
   "display_name": "Python 3.10.2 ('geo3')",
   "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.10.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}