Deployed 39b946f with MkDocs version: 1.5.3

03_compile_time_series/03_compile_time_series.ipynb

@@ -0,0 +1,228 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# **Compile Time Series by Units of Analysis**\n",
+    "\n",
+    "Authors: \n",
+    "- Victor Tang\n",
+    "\n",
+    "Reviewed/Edited by:\n",
+    "- Marcos Kavlin\n",
+    "- Dr. Andrew Dean\n",
+    "\n",
+    "### Purpose\n",
+    "\n",
+    "This notebook is the third in the Wetland Function Assessment Tutorial.\n",
+    "The goal of this notebook is to show you, the user, how to aggregate your the time series information by the landscape units of analysis that were calculated in the previous notebook.\n",
+    "\n",
+    "\n",
+    "### Workflow\n",
+    "\n",
+    "1. Import required packages.\n",
+    "2. Load Sentinel 1 imagery.\n",
+    "3. Load landscape unit polygons.\n",
+    "4. Aggregate the Sentinel 1 time series pixel information by landscape unit.\n",
+    "5. Export results.\n",
+    "\n",
+    "### Notes\n",
+    "\n",
+    "The composites that were used to run this notebook are placeholders used to demonstrate this workflow. The data required to run the steps demonstrated in this notebook are 10 day Sentinel 1 median composites. In this tutorial the images were produced for the year 2022. The images were grouped into 10 day composites as that was most conducive to a complete year-long time series, while still removing speckle and noise. They were all stored in a folder for the year 2022.\n",
+    "\n",
+    "The landscape units used in this tutorial were the output from the previous notebook."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Note Description:**\n",
+    "\n",
+    "This notebook aggregate pixels in Sentinel-1 image time series by landscape units, which were extracted from segmentation of Sentinel-2 composite. It returns an \"average\" time series for each of landscape units that represent the general temporal pattern of S1 backscatter for a certain year. The temporal resolution of the \"average\" time sereis is determined by the interval of Sentinel-1 images, which is 10 days in this case."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 1. Import required packages"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "import glob\n",
+    "import geopandas as gpd\n",
+    "\n",
+    "import xarray as xr\n",
+    "import rioxarray as rxr\n",
+    "\n",
+    "from rasterio.features import geometry_mask\n",
+    "\n",
+    "base_dir = f\"{path_to_data}/\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 2. Load Sentinel-1 data\n",
+    "\n",
+    "The code below was written with the imagery named using the following naming convention: \"*an_giang_2022_10d_**02**.tif*\"\n",
+    "  - In this example 'an_giang' is the province name\n",
+    "  - 2022 is the year in question\n",
+    "  - 10d specifies the timespan of the composite, 02 specifies what it's order (temporally) is within the year.\n",
+    "\n",
+    "The code completes the following steps:\n",
+    "  - Lists all the .tif files in the specified directory\n",
+    "  - Make a list of the order numbers of each image\n",
+    "  - Open each image, select the VH band and append it to an empty list\n",
+    "  - Concatenate da_list, using the list of order numbers as a pandas.Index, for the dimension argument.\n",
+    "  - Convert the units of Sentinel-1 values from linear to decibel units (for better visualization)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 84,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# find all the .tif files in a given folder\n",
+    "fpath_list = glob.glob(base_dir + \"*.tif\")\n",
+    "\n",
+    "# extract image order number from file name\n",
+    "x = [int(item.split(\"_\")[-1].split(\".\")[0]) for item in fpath_list]\n",
+    "\n",
+    "# load all the tif files as a list of xarray DataArray\n",
+    "da_list = []\n",
+    "for fpath in fpath_list:\n",
+    "    da = rxr.open_rasterio(fpath)  # load Sentinel-1 image\n",
+    "    da = da.sel(band=1)  # select VH band\n",
+    "    da_list.append(da)\n",
+    "\n",
+    "# concatenate the list of xarray DataArray\n",
+    "da = xr.concat(da_list, pd.Index(x))\n",
+    "\n",
+    "# convert unit of Sentinel-1 values from linear to decibel\n",
+    "da.values = 10 * np.log10(da.values)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 3. Load polygons of landscape units\n",
+    "\n",
+    "Once we have loaded our Sentinel-1 data the next step is to load the landscapes units we need to aggregate the data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 85,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf = gpd.read_file(\"data/an_giang_segmentation.geojson\")\n",
+    "gdf = gdf.to_crs(da.rio.crs.to_epsg())\n",
+    "gdf = gdf.set_index(\"PID\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 4. Aggregate pixels by landscape units\n",
+    "\n",
+    "In order to aggregate our pixels by landscape unit, we need to iterate through the features in our landscape unit geodataframe.\n",
+    "\n",
+    "The code below executes the following steps:\n",
+    "- Obtain the number of time steps in our time series.\n",
+    "- Create an empty pandas.DataFrame, in which to store the aggregated values for each landscape unit.\n",
+    "- Iterate through the landscape units:\n",
+    "  - Clip the time series data, by the spatial extent of thelandscape unit.\n",
+    "  - Obtain the median of the remaining pixels\n",
+    "  - Generate a 1-dimensional array of medians, which will be used as the time series for each landscape unit.\n",
+    "  - Append this time series to the appropriate index in our empty pandas.DataFrame."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 86,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n = da.shape[0]  # number of images (i.e. time steps)\n",
+    "df = pd.DataFrame(index=gdf.index, columns=np.arange(1, n + 1))\n",
+    "\n",
+    "for pid, row in gdf.iterrows():\n",
+    "    # subset image to improve efficiency\n",
+    "    xmin, ymin, xmax, ymax = row.geometry.bounds\n",
+    "    da2 = da.sel(x=slice(xmin, xmax), y=slice(ymax, ymin))\n",
+    "\n",
+    "    mask = geometry_mask(\n",
+    "        gpd.GeoSeries([row.geometry]),\n",
+    "        out_shape=da2.shape[-2:],  # get dimension of x and y\n",
+    "        transform=da2.rio.transform(),\n",
+    "    )\n",
+    "    mask = ~mask\n",
+    "\n",
+    "    data = da2.values[:, mask]  # select pixels within polygon mask\n",
+    "    df.loc[pid] = np.median(data, axis=-1)  # aggregate by median"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 5. Export the resulting DataFrame\n",
+    "\n",
+    "Once we've filled our new DataFrame with a median time series for each landscape unit, we save it, to use in the next notebook.\n",
+    "\n",
+    "***This ouput data is provided so once can run the following notebooks***"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df.to_csv(f\"{output_path}/vh_2022.csv\")"
+   ]
+  }
+ ],
+ "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.10.11"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "80ea545c89f4740f2190761fe3f09035380c8ff11000f9cb62822d0fef0270af"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}