MODIS_workflow_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "938a2879-f85c-4fc6-bb28-f91d6a680644",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    },
    "tags": []
   },
   "source": [
    "### This notebook contains code for the processing of raw MODIS geotiff files as used and described in Lang et al. (xxxx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ec51085-bab2-47e9-bdf3-d3030bb4b656",
   "metadata": {},
   "source": [
    "Inputs: \n",
    "1) MODIS fSCA (from MODSCAG) and RF (from MODDRFS) daily.tif files in default sinusoidal projection.\n",
    "2) Shapefile for watershed/region of interest.\n",
    "\n",
    "Outputs:\n",
    "1) Step - Folder with MODIS data clipped by watershed.\n",
    "2) Step - Folder with MODIS data clipped by watershed and reprojected.\n",
    "3) Step - Folder with binary snow cover fraction mask rasters.\n",
    "4) Step - Dataset containing data from clipped and reprojected RF, fSCA, and mask rasters.\n",
    "5) Step - Same dataset as Step 4, but also containing RF masked by snow cover fraction.\n",
    "6) Step (final output) - Same dataset as Step 5, but containing RF masked by snow cover fraction and snow persistence.\n",
    "\n",
    "Note: steps 1-2 must be completed for both RF and fSCA data seperately. Steps 4-6 must be completed in order and rely on outputs from steps 1-3.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d85966ba-9573-4667-8b0f-00a6427767ff",
   "metadata": {},
   "source": [
    "###### Author: Otto Lang. Assited by: Patrick Naple and Dillon Ragar.\n",
    "###### Contact otto.lang@utah.edu for questions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18a7a278-dc45-4bcb-b980-36d4c5dca871",
   "metadata": {},
   "source": [
    "#### 1. Clipping "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f73208a0-f622-41f4-8731-a43c8f359851",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Clipping daily MODIS Geotiffs to the watershed of interest \n",
    "#Creates daily clipped geotiffs and stores them in an output directory\n",
    "\n",
    "#Loading packages for this step\n",
    "import gdal\n",
    "import fnmatch\n",
    "import os\n",
    "gdal.UseExceptions()\n",
    "\n",
    "#Directories (update this)\n",
    "inDir = \"/Path/to/folder/with/MODIS/geotiffs\"\n",
    "outDir = \"/Path/to/clipped/output/folder\"\n",
    "shpin = \"/Path/to/watershed/shapefile.shp\"\n",
    "os.chdir (inDir)\n",
    "\n",
    "#Defining findRasters function\n",
    "def findRasters (path, filter):\n",
    "    for root, dirs, files in os.walk(path, filter):\n",
    "        for file in fnmatch.filter(files, filter):\n",
    "            yield os.path.join (root, file)\n",
    "\n",
    "#Creating new rasters in outDir directory\n",
    "for raster in findRasters (inDir, '*.tif'):\n",
    "    (infilepath, infilename)= os.path.split (raster)\n",
    "    print(infilename)\n",
    "    outraster = outDir+ 'clip_'+ infilename\n",
    "    print(outraster)\n",
    "    result = gdal.Warp(outraster, raster, cutlineDSName = shpin, cropToCutline=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "963026a5-34e8-48ed-9629-78e6e2f248d9",
   "metadata": {},
   "source": [
    "#### 2. Reprojecting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3bc2b03-ba2f-4ff7-813a-d2efae75d33d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reprojecting clipped daily MODIS Geotiffs from default sinusoidal projection to desired projection (WGS 1984 UTM 12N shown)\n",
    "# This code creates a new directory with reprojected clipped geotiffs\n",
    "\n",
    "#Loading packages for this step\n",
    "import gdal\n",
    "import fnmatch\n",
    "import os\n",
    "gdal.UseExceptions()\n",
    "\n",
    "#Directories (update this)\n",
    "inDir = \"/Path/to/folder/with/clipped/geotiffs\"\n",
    "outDir = \"/Path/to/reprojected/output/folder\"\n",
    "os.chdir (inDir)\n",
    "\n",
    "#Defining findRasters function (no need if done in previous step)\n",
    "def findRasters (path, filter):\n",
    "    for root, dirs, files in os.walk(path, filter):\n",
    "        for file in fnmatch.filter(files, filter):\n",
    "            yield os.path.join (root, file)\n",
    "            \n",
    "#Creating new rasters in outDir directory\n",
    "for raster in findRasters (inDir, '*.tif'):\n",
    "    (infilepath, infilename)= os.path.split (raster)\n",
    "    print(infilename)\n",
    "    outraster = outDir+ 'reprj_'+ infilename\n",
    "    print(outraster)\n",
    "    # Edit below for desired projection\n",
    "    result = gdal.Warp(outraster, raster, srcSRS = '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs',\n",
    "                   dstSRS = '+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d8b3020-507b-4e45-9f96-3f960711fe29",
   "metadata": {},
   "source": [
    "#### 3. Creating binary mask rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9250c461-94ba-4511-8842-9716f9f73693",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creating snow cover fraction mask rasters for each day\n",
    "#Creates daily binary rasters where values greater than the snow cover fraction threshold = 1, and others = 0. \n",
    "\n",
    "#Loading packages for this step\n",
    "import gdal\n",
    "import re\n",
    "import os\n",
    "import fnmatch\n",
    "from matplotlib import pyplot\n",
    "import numpy as np\n",
    "\n",
    "#Directories (update this)\n",
    "inDirscf = \"/Path/to/folder/with/clipped/reprojected/snow/cover/geotiffs\"\n",
    "outDirscf = \"/Path/to/output/snow/cover/mask/folder/\"\n",
    "\n",
    "#Defining functions\n",
    "def findRasters (path, filter):\n",
    "    for root, dirs, files in os.walk(path, filter):\n",
    "        for file in fnmatch.filter(files, filter):\n",
    "            yield os.path.join (root, file)\n",
    "\n",
    "def CreateGeoTiff(outRaster, data, geo_transform, projection):\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    rows, cols = data.shape\n",
    "    DataSet = driver.Create(outRaster, cols, rows, 1, gdal.GDT_Byte)\n",
    "    DataSet.SetGeoTransform(geo_transform)\n",
    "    DataSet.SetProjection(projection)\n",
    "    DataSet.GetRasterBand(1).WriteArray(data)\n",
    "    DataSet.GetRasterBand(1).SetNoDataValue(0)\n",
    "    DataSet = None\n",
    "    \n",
    "#Creating mask rasters   \n",
    "for raster in findRasters (inDirscf, '*.tif'):\n",
    "    (infilepath, infilename)= os.path.split (raster)\n",
    "    print(infilename)\n",
    "    scf = gdal.Open(raster, gdal.GA_ReadOnly)\n",
    "    gt = scf.GetGeoTransform()\n",
    "    proj = scf.GetProjection()\n",
    "    scf_band = scf.GetRasterBand(1)\n",
    "    scf_data = scf_band.ReadAsArray()\n",
    "    #Edit fSCA values for desired snow cover fraction mask\n",
    "    scf_mask = np.where(np.logical_and(scf_data>=95, scf_data<=100)==True, 1, 0)\n",
    "    scf_mask_raster = outDirscf+ infilename\n",
    "    data = scf_mask\n",
    "    geo_transform = gt\n",
    "    projection = proj\n",
    "    CreateGeoTiff(scf_mask_raster, data, geo_transform, projection)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4c574c7-73c0-49b0-9119-dbbb91511ce8",
   "metadata": {},
   "source": [
    "#### 4. Creating dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ea24cb4-4e26-4a75-b304-884d97f48c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Creating xarray dataset for all clipped and reprojected RF, fSCA, and mask rasters \n",
    "\n",
    "#Loading packages for this step\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import dask\n",
    "from datetime import datetime\n",
    "import rioxarray\n",
    "\n",
    "#Directories (update this)\n",
    "mask_files = \"/Path/to/folder/with/mask/rasters/*.tif\"\n",
    "scf_values = \"/Path/to/folder/with/clipped/and/reprojected/snow/cover/rasters/*.tif\"\n",
    "rf_values = \"/Path/to/folder/with/clipped/and/reprojected/RF/rasters/*.tif\"\n",
    "\n",
    "#Defining function to add time dimension to dataset\n",
    "def add_time_dim(xda):\n",
    "    xda = xda.expand_dims(time = [datetime.now()])\n",
    "    return xda\n",
    "\n",
    "#Adding RF data \n",
    "rf = xr.open_mfdataset(rf_values, preprocess = add_time_dim)\n",
    "\n",
    "#Adding dates (Update for dates of timeseries)\n",
    "dates = pd.date_range('2000-02-01', '2023-09-30', freq='D')\n",
    "yrs = []\n",
    "yrs.append(dates)                          \n",
    "t_range = yrs[0].union_many(yrs[1:]) \n",
    "rf['time'] = t_range\n",
    "\n",
    "#Adding mask data \n",
    "mask = xr.open_mfdataset(mask_files, preprocess = add_time_dim)\n",
    "mask['time'] = t_range\n",
    "\n",
    "#Adding fSCA data \n",
    "scf_val = xr.open_mfdataset(mask_values, preprocess = add_time_dim)\n",
    "scf_val['time'] = t_range\n",
    "\n",
    "#Assigning variables to dataset \n",
    "data = rf.assign(rf = rf['band_data'], mask=mask['band_data'], scf=scf_val['band_data'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3021302e-b281-45f2-a1ab-7c5324051b94",
   "metadata": {},
   "source": [
    "#### 5. Snow cover fraction mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69cf3393-20ce-4e4b-a0fd-cdaf34ebb985",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Snow cover fraction mask step\n",
    "\n",
    "#Adding 1 to retain RF values of 0 through mask calculation step\n",
    "data['rf'] = data['rf']+1\n",
    "#Mask calculation\n",
    "data['rf_masked'] = data['rf'] * data['mask']\n",
    "#Removing masked out data\n",
    "data['rf_masked'] = data['rf_masked'].where(data['rf_masked'] != 0)\n",
    "#Subtraction 1 to revert back to true RF values\n",
    "data['rf_masked'] = data['rf_masked']-1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c850aba5-f081-49e5-8b9c-29f810536306",
   "metadata": {},
   "source": [
    "#### 6. Snow persistence mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e80c67a-d720-44a9-8ad2-c8d5c6ad9eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Snow persistence mask step\n",
    "\n",
    "#Selecting only April-June (This is the \"snowmelt season\" - see paper)\n",
    "def is_amj(month):\n",
    "    return (month >= 4) & (month <= 6)\n",
    "data['rf_daily_AMJ'] = data['rf_masked'].sel(time=is_amj(data['time.month']))\n",
    "\n",
    "#Creating pixelcount variable: if there is data in a pixel- replace with 1, if not = replace with 0\n",
    "data['pixelcount'] = data['rf_daily_AMJ'].where(data['rf_daily_AMJ']>=0, 0)\n",
    "data['pixelcount'] = data['pixelcount'].where(data['pixelcount']==0, 1)\n",
    "\n",
    "#Creating variable showing average percent (from 0 to 1) that a given pixel is included over the timeframe\n",
    "data['average_pixelcount'] = data['pixelcount'].mean(dim = 'time', skipna = True)\n",
    "\n",
    "#Creating single snow persistence mask where snow persistence >= 5% - replaced with 1, otherwise - replaced with 0\n",
    "data['pixelmask'] = data['average_pixelcount'].where(data['average_pixelcount']>=.05, 0)\n",
    "data['pixelmask'] = data['pixelmask'].where(data['pixelmask']<.05, 1)\n",
    "\n",
    "#Multiplying RF data by snow persistence mask (making sure to include 0s by adding 1 to the data, and removing 1 again after the masking)\n",
    "data['rf_daily_AMJ'] = data['rf_daily_AMJ']+1\n",
    "data['rf_final'] = data['pixelmask'] * data['rf_daily_AMJ']\n",
    "data['rf_final'] = data['rf_final'].where(data['rf_final'] != 0)\n",
    "data['rf_final'] = data['rf_final']-1\n",
    "\n",
    "#data['rf_final'] is the fully processed RF data quantified in the paper. "
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}