PS_preprocess_reflectance_images.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "52ce4a11",
   "metadata": {},
   "source": [
    "### Planetscope preprocessing - Reflectance Images"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46a9410c",
   "metadata": {},
   "source": [
    "This script:\n",
    "- Reads multiple Planetscope images/scenes   \n",
    "\n",
    "- Convert TOA Radiance to TOA Reflectance   \n",
    "\n",
    "- Clips images over selected area of interest (AOI).   \n",
    "\n",
    "Main packages used are rasterio and fiona.  \n",
    "\n",
    "Code is based on notebooks from Planetlabs, available here: https://github.com/planetlabs/notebooks  \n",
    "\n",
    "Data used are from Planetscope 21.-26. June 2019 over a selected study area in the Brazilian Amazon.   \n",
    "\n",
    "Available export options:\n",
    "- Reflectance scaled - *(default setting)*\n",
    "- Reflectance scaled clipped to AOI - *(default setting)*   \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "- Reflectance unscaled\n",
    "- Reflectance unscaled clipped to AOI\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a6407d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install libraries\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import rasterio\n",
    "import fiona"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b470d4ea",
   "metadata": {},
   "source": [
    "#### Step 1. Set up project folders - rest is automatic\n",
    "The following input and directory structure is required:\n",
    "\n",
    "- *Number of images/scenes*\n",
    "\n",
    "\n",
    "- *Project folder*  -  all folders inside Project folder\n",
    "\n",
    "\n",
    "- *AnalyticMS folder*  -  only the AnalyticMS.tif files, no unusable datamasks(udm)\n",
    "\n",
    "\n",
    "- *Metadata folder*  -  only the metadata.xml files\n",
    "\n",
    "\n",
    "- *AOI folder*  -  just one AOI file in GeoJSON format\n",
    "\n",
    "\n",
    "- *Results folder*  -  receives export of reflectance images, also serves as input to the clip-functions that further clips the results to selected AOI\n",
    "\n",
    "\n",
    "The default setting export scaled images and scaled images clipped to AOI.\n",
    "\n",
    "The unscaled images and unscaled clipped images are available by calling the *unscaled-functions* in Step 6, then Step 8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "800f81f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_images = 6\n",
    "Project_folder = 'planet_21-26/'\n",
    "AnalyticMS_folder = 'planet_21-26/AnalyticMS_files/'\n",
    "Metadata_folder = 'planet_21-26/Metadata_files/'\n",
    "aoi_folder = 'planet_21-26/AOI/'\n",
    "Results_folder = 'planet_21-26/Results/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "865d5d88",
   "metadata": {},
   "source": [
    "#### Step 2. Read and Inspect Images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbe3632c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "planet_21-26/AnalyticMS_files/20190721_133242_1034_3B_AnalyticMS.tif\n",
      "planet_21-26/AnalyticMS_files/20190722_133028_103a_3B_AnalyticMS.tif\n",
      "planet_21-26/AnalyticMS_files/20190723_131552_0f3d_3B_AnalyticMS.tif\n",
      "planet_21-26/AnalyticMS_files/20190724_131505_1_1050_3B_AnalyticMS.tif\n",
      "planet_21-26/AnalyticMS_files/20190725_133328_1042_3B_AnalyticMS.tif\n",
      "planet_21-26/AnalyticMS_files/20190726_133225_103a_3B_AnalyticMS.tif\n"
     ]
    }
   ],
   "source": [
    "# Use 6 scenes from 21.-26. June 2019 for coversion to TOA Reflectance\n",
    "# Collect image filenames \n",
    "from os import walk\n",
    "AnalyticMS_files = []\n",
    "for (dirpath, dirnames, filenames) in walk(AnalyticMS_folder):\n",
    "    AnalyticMS_files.extend(filenames)\n",
    "    break\n",
    "\n",
    "# Add image file path\n",
    "AnalyticMS_path = []\n",
    "for i in AnalyticMS_files:\n",
    "    AnalyticMS_path.append(AnalyticMS_folder + i)\n",
    "    print(AnalyticMS_folder + i) # Verify correct images"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6461646",
   "metadata": {},
   "source": [
    "Note the ordering of the images: 1,2,3,4,5,6 for dates: 21,22,23,24,25,26"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "505a2e22",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "planet_21-26/Metadata_files/20190721_133242_1034_3B_AnalyticMS_metadata.xml\n",
      "planet_21-26/Metadata_files/20190722_133028_103a_3B_AnalyticMS_metadata.xml\n",
      "planet_21-26/Metadata_files/20190723_131552_0f3d_3B_AnalyticMS_metadata.xml\n",
      "planet_21-26/Metadata_files/20190724_131505_1_1050_3B_AnalyticMS_metadata.xml\n",
      "planet_21-26/Metadata_files/20190725_133328_1042_3B_AnalyticMS_metadata.xml\n",
      "planet_21-26/Metadata_files/20190726_133225_103a_3B_AnalyticMS_metadata.xml\n"
     ]
    }
   ],
   "source": [
    "# Reading image metadata\n",
    "Analytic_meta_files = []\n",
    "for (dirpath, dirnames, filenames) in walk(Metadata_folder):\n",
    "    Analytic_meta_files.extend(filenames)\n",
    "    break\n",
    "\n",
    "# Add image metadata file path\n",
    "Analytic_meta_path = []\n",
    "for i in Analytic_meta_files:\n",
    "    Analytic_meta_path.append(Metadata_folder + i)\n",
    "    print(Metadata_folder + i) # Verify correct metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "61b6cac3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of img_list: 6 \n"
     ]
    }
   ],
   "source": [
    "# Loading images into list using rasterio\n",
    "img_list = []\n",
    "for image in AnalyticMS_path:\n",
    "    img_list.append(rasterio.open(image))\n",
    "print('Length of img_list: {} '.format(len(img_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9375c0c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image: dtype | crs | band count\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n"
     ]
    }
   ],
   "source": [
    "# Inspecting image metadata\n",
    "print(\"Image: dtype | crs | band count\")\n",
    "for image in img_list:\n",
    "    print(image.meta['dtype'], image.meta['crs'], image.meta['count'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7ee7c2d",
   "metadata": {},
   "source": [
    "Let's take a closer look at what these bands contain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "68c13e7a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blue\n",
      "green\n",
      "red\n",
      "undefined\n"
     ]
    }
   ],
   "source": [
    "# Read in color interpretations of each band in img1 - assume same values for rest of the images\n",
    "colors = [img_list[0].colorinterp[band] for band in range(img_list[0].count)]\n",
    "\n",
    "# taking a look at img1's band types:\n",
    "for color in colors:\n",
    "    print(color.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62f9df3e",
   "metadata": {},
   "source": [
    "The fourth channel is a binary alpha mask: this is common in satellite color models, and can be confirmed in Planet's documentation on the PSSCene3Band product."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e19a816d",
   "metadata": {},
   "source": [
    "#### Step 3. Extract the Data from Each Spectral Band\n",
    "In this step, Rasterio (a Python library for reading and writing geospatial raster datasets) is used to open the raster images (the .tif files). \n",
    "\n",
    "Then, the band data will is extracted and loaded into arrays for further manipulation with Python's NumPy libary.\n",
    "\n",
    "Note: in PlanetScope 4-band images, the band order is BGRN: (1) Blue, (2) Green, (3) Red, (4) Near-infrared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "76e4ace6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Radiance values are loaded into lists\n",
    "band_blue_radiance = [0]*n_images\n",
    "band_green_radiance = [0]*n_images\n",
    "band_red_radiance = [0]*n_images\n",
    "band_nir_radiance = [0]*n_images\n",
    "\n",
    "# Using rasterio to read radiance values for the images\n",
    "for i in range(n_images):\n",
    "    with rasterio.open(AnalyticMS_path[i]) as src:\n",
    "        band_blue_radiance[i] = src.read(1)\n",
    "    with rasterio.open(AnalyticMS_path[i]) as src:\n",
    "        band_green_radiance[i] = src.read(2)\n",
    "    with rasterio.open(AnalyticMS_path[i]) as src:\n",
    "        band_red_radiance[i] = src.read(3)\n",
    "    with rasterio.open(AnalyticMS_path[i]) as src:\n",
    "        band_nir_radiance[i] = src.read(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35aac229",
   "metadata": {},
   "source": [
    "#### Step 4. Extract the Coefficients\n",
    "Before converting to reflectance, the conversion coefficients from the metadata file (the .xml file) must be extracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1ec7c24c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1: 2.22957954184e-05, 2: 2.36158292355e-05, 3: 2.63003228576e-05, 4: 3.95838798777e-05}\n",
      "{1: 2.22520864862e-05, 2: 2.35908050459e-05, 3: 2.63407248495e-05, 4: 3.9577582354e-05}\n",
      "{1: 2.3376622001e-05, 2: 2.47621088261e-05, 3: 2.76131113868e-05, 4: 4.14097363863e-05}\n",
      "{1: 2.3524916513e-05, 2: 2.48470408363e-05, 3: 2.77309803512e-05, 4: 4.13036148082e-05}\n",
      "{1: 2.21525543854e-05, 2: 2.34361488133e-05, 3: 2.61307246766e-05, 4: 3.90849967045e-05}\n",
      "{1: 2.20519050955e-05, 2: 2.33785804456e-05, 3: 2.61037617704e-05, 4: 3.92215395407e-05}\n"
     ]
    }
   ],
   "source": [
    "from xml.dom import minidom\n",
    "\n",
    "# Gathering the coefficients for the 7 images\n",
    "coeffs_list = [] \n",
    "for image in Analytic_meta_path:\n",
    "    xmldoc = minidom.parse(image)\n",
    "    nodes = xmldoc.getElementsByTagName(\"ps:bandSpecificMetadata\")\n",
    "    \n",
    "    # XML parser refers to bands by numbers 1-4\n",
    "    coeffs = {} # Coefficients for each image/scene\n",
    "    for node in nodes:\n",
    "        band_nr = node.getElementsByTagName(\"ps:bandNumber\")[0].firstChild.data\n",
    "        if band_nr in ['1', '2', '3', '4']:\n",
    "            i = int(band_nr)\n",
    "            value = node.getElementsByTagName(\"ps:reflectanceCoefficient\")[0].firstChild.data\n",
    "            coeffs[i] = float(value)\n",
    "    coeffs_list.append(coeffs)\n",
    "\n",
    "for coffee in coeffs_list:\n",
    "    print (coffee) # Verify data/coffee"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07604cbe",
   "metadata": {},
   "source": [
    "Note that the coefficients are all of order 1e<sup>-5</sup>, and that the coefficient for NIR is significantly higher than the coefficient for blue. This is a big deal if your use case involves performing band math because a pixel with a NIR/blue ratio of 1.0 in the radiance image will have a NIR/blue ratio of 3.35/1.929=1.73 in the reflectance image.   \n",
    "Most spectral indices are defined in terms of reflectance, not radiance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77a9178e",
   "metadata": {},
   "source": [
    "#### Step 5: Convert Radiance to Reflectance\n",
    "Radiance is measured in SI units: $W/m^2$. Reflectance is a ratio from 0 to 1. The conversion is performed as a per-band scalar multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6325f9b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red band radiance goes from: 0 to 12546\n",
      "Red band reflectance goes from: 0.0 to 0.3299638505714496\n"
     ]
    }
   ],
   "source": [
    "# Radiance values are loaded into lists\n",
    "band_blue_reflectance = [];   blue_ref_scaled = [0]*n_images\n",
    "band_green_reflectance = [];  green_ref_scaled = [0]*n_images\n",
    "band_red_reflectance = [];    red_ref_scaled = [0]*n_images\n",
    "band_nir_reflectance = [];    nir_ref_scaled = [0]*n_images\n",
    "\n",
    "# Calculating reflectance from radiance and conversion coefficients\n",
    "for i in range(n_images):\n",
    "        band_blue_reflectance.append(band_blue_radiance[i] * coeffs_list[i][1])\n",
    "        band_green_reflectance.append(band_green_radiance[i] * coeffs_list[i][2])\n",
    "        band_red_reflectance.append(band_red_radiance[i] * coeffs_list[i][3])\n",
    "        band_nir_reflectance.append(band_nir_radiance[i] * coeffs_list[i][4])\n",
    "        \n",
    "# Check the results by looking at the min-max range of the radiance (before) vs reflectance (after)\n",
    "# for the red band from the first image, excluding NaN values\n",
    "red_rad_min = np.nanmin(band_red_radiance[0]);    red_rad_max = np.nanmax(band_red_radiance[0])\n",
    "red_ref_min = np.nanmin(band_red_reflectance[0]); red_ref_max = np.nanmax(band_red_reflectance[0])\n",
    "print(\"Red band radiance goes from: {} to {}\".format(red_rad_min, red_rad_max))\n",
    "print(\"Red band reflectance goes from: {} to {}\".format(red_ref_min, red_ref_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c0e22f9",
   "metadata": {},
   "source": [
    "#### Step 6. Write the Reflectance Image\n",
    "A note: Reflectance is generally defined as a floating point number between 0 and 1, but image file formats are much more commonly stored as unsigned integers. A common practice in the industry is to multiply the reflectance value by a *scale factor* of 10,000, then save the result as a file with data type uint16.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa35302e",
   "metadata": {},
   "source": [
    "Export of normal scaled GeoTIFF (multiplied with scale factor of 10 000 and stored as *dtype=uint16*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "74f02f10",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to export scaled images - (alternative option)\n",
    "def scaled_export(images):\n",
    "    for i in range(len(images)):\n",
    "        # Get the metadata of original GeoTIFF:\n",
    "        meta = images[i].meta\n",
    "        \n",
    "        # Set the source metadata as kwargs we'll use to write the new data:\n",
    "        # update the 'dtype' value to uint16:\n",
    "        kwargs = meta\n",
    "        kwargs.update(dtype='uint16')\n",
    "        \n",
    "        # As noted above, scale reflectance value by a factor of 10k:\n",
    "        scale = 10000\n",
    "        blue_ref_scaled[i] = scale * band_blue_reflectance[i]\n",
    "        green_ref_scaled[i] = scale * band_green_reflectance[i]\n",
    "        red_ref_scaled[i] = scale * band_red_reflectance[i]\n",
    "        nir_ref_scaled[i] = scale * band_nir_reflectance[i]\n",
    "        \n",
    "        # Compute new min & max values for the scaled red band, just for comparison\n",
    "        red_min_scaled = np.amin(red_ref_scaled[0])\n",
    "        red_max_scaled = np.amax(red_ref_scaled[0])\n",
    "        \n",
    "        # Convert to type 'uint16'\n",
    "        from rasterio import uint16\n",
    "        blue_scaled = blue_ref_scaled[i].astype(uint16)\n",
    "        green_scaled = green_ref_scaled[i].astype(uint16)\n",
    "        red_scaled = red_ref_scaled[i].astype(uint16)\n",
    "        nir_scaled = nir_ref_scaled[i].astype(uint16)\n",
    "        \n",
    "        # New name for exported image\n",
    "        img_name_ref_scaled = (Results_folder + AnalyticMS_files[i]).replace('.tif','_Ref.tif')\n",
    "        \n",
    "        # Write band calculations to a new unscaled GeoTIFF file with same band order (BGRN)\n",
    "        with rasterio.open(img_name_ref_scaled, 'w', **kwargs) as dst:\n",
    "                dst.write_band(1, blue_scaled)\n",
    "                dst.write_band(2, green_scaled)\n",
    "                dst.write_band(3, red_scaled)\n",
    "                dst.write_band(4, nir_scaled)\n",
    "         # Comparing min & max values before/after scaling, with the first image red band\n",
    "        if i == 0:\n",
    "            print(\"With scale factor: {}\".format(scale))\n",
    "            print(\"Before scaling: Red band reflectance from: {} to {}\"\\\n",
    "                  .format(red_ref_min, red_ref_max))\n",
    "            print(\"After scaling: Red band reflectance from: {} to {}\"\\\n",
    "                  .format(red_min_scaled,red_max_scaled))\n",
    "        if i == (len(images)-1):\n",
    "            print(\"Success, {} scaled images(dtype={}) have been exported to your Results folder: {}\"\\\n",
    "                  .format(n_images, kwargs['dtype'], Results_folder ))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "926e017e",
   "metadata": {},
   "source": [
    "Export of unscaled GeoTIFF, dtype=float64 with Reflectance values between 0.0 - 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ad09d10e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to export unscaled images - (alternative option)\n",
    "def unscaled_export(images):\n",
    "    for i in range(len(images)):\n",
    "        # Get the metadata of original GeoTIFF:\n",
    "        meta = images[i].meta\n",
    "\n",
    "        # Set the source metadata as kwargs we'll use to write the new data:\n",
    "        # update the 'dtype' value to float64:\n",
    "        kwargs = meta\n",
    "        kwargs.update(dtype='float64')\n",
    "\n",
    "        # New name for exported image\n",
    "        img_name_ref_unscaled = (Results_folder + AnalyticMS_files[i]).replace('.tif','_Ref_unscaled.tif')\n",
    "\n",
    "        # Write band calculations to a new unscaled GeoTIFF file with same band order (BGRN)\n",
    "        with rasterio.open(img_name_ref_unscaled, 'w', **kwargs) as dst:\n",
    "                dst.write_band(1, band_blue_reflectance[i])\n",
    "                dst.write_band(2, band_green_reflectance[i])\n",
    "                dst.write_band(3, band_red_reflectance[i])\n",
    "                dst.write_band(4, band_nir_reflectance[i])\n",
    "        \n",
    "        if i == (len(images)-1):\n",
    "            print(\"Success, {} unscaled images(dtype={}) have been exported to your Results folder: {}\"\\\n",
    "                  .format(n_images, kwargs['dtype'], Results_folder ))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ed666a8",
   "metadata": {},
   "source": [
    "##### Export Reflectance Images to your Results folder\n",
    "- For the scaled images, run the function *scaled_export()* with input(*img_list*)   \n",
    "\n",
    "\n",
    "- For the unscaled images, run the function *unscaled_export()* with input(*img_list*)\n",
    "\n",
    "Note that unscaled(*dtype=float64*) takes 4 times as much storage space as scaled(*dtype=uint16*)\n",
    "\n",
    "For clipping reflectance to AOI, proceed further"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "47d7fb8b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With scale factor: 10000\n",
      "Before scaling: Red band reflectance from: 0.0 to 0.3299638505714496\n",
      "After scaling: Red band reflectance from: 0.0 to 3299.638505714496\n",
      "Success, 6 scaled images(dtype=uint16) have been exported to your Results folder: planet_21-26/Results/\n"
     ]
    }
   ],
   "source": [
    "# Uncomment and run for wanted image type:\n",
    "# Export scaled images:\n",
    "scaled_export(img_list)\n",
    "\n",
    "# Export unscaled images:\n",
    "#unscaled_export(img_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e7cee5ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load image results into lists for the clip functions\n",
    "Results_files = []\n",
    "Results_path_scaled = []; Results_path_unscaled = []\n",
    "img_list_scaled = []; img_list_unscaled = []\n",
    "\n",
    "for (dirpath, dirnames, filenames) in walk(Results_folder):\n",
    "    Results_files.extend(filenames)\n",
    "    break\n",
    "\n",
    "for filename in Results_files:\n",
    "    if \"Ref.tif\" in filename:\n",
    "        Results_path_scaled.append(Results_folder + filename)\n",
    "        img_list_scaled.append(rasterio.open(Results_folder + filename))\n",
    "    if \"Ref_unscaled.tif\" in filename:\n",
    "        Results_path_unscaled.append(Results_folder + filename)\n",
    "        img_list_unscaled.append(rasterio.open(Results_folder + filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e310f619",
   "metadata": {},
   "source": [
    "Using rasterio.plot (a matplotlib interface) to preview the results of our mosaic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "311ed615",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show\n",
    "#show(img_list[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dfc6500",
   "metadata": {},
   "source": [
    "#### Step 7. Clip Images to AOI Boundaries\n",
    "A mask is used to clip the images to the extents of our AOI.\n",
    "\n",
    "For this operation rasterio's sister-library fiona will read in our AOI (as a GeoJSON file).\n",
    "Just as rasterio is used to manipulate raster data, fiona works similarly on vector data. Where rasterio represents raster imagery as numpy arrays, fiona represents vector data as GeoJSON-like Python dicts.\n",
    "\n",
    "After the GeoJSON is read, the geometry of the AOI is extracted, with the geometry as dict key."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19836261",
   "metadata": {},
   "source": [
    "##### A note about Coordinate Reference Systems\n",
    "\n",
    "Before the AOI can be applied to the images, the Coordinate Reference System (CRS) need to match.\n",
    "This can be checked by reading the crs attribute of the Collection object generated by *fiona.open()*.\n",
    "\n",
    "In this example, the CRS of the AOI is: *EPSG:4326*.\n",
    "\n",
    "While the CRS for the images is: *EPSG:32722*.\n",
    "\n",
    "\n",
    "Before the clip can be applied, the geometry of the AOI needs to be transformed to match the CRS of the imagery.\n",
    "Luckily, fiona is smart enough to apply the necessary mathematical transformation to a set of coordinates in order to convert them to new values: \n",
    "apply *fiona.transform.transform_geom* to the AOI geometry to do this, specifying the GeoJSON's CRS as the source CRS, and the imagery's CRS as the destination CRS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "49f46b87",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "planet_21-26/AOI/AOI_small.geojson\n"
     ]
    }
   ],
   "source": [
    "# Reading AOI filename - (assumes just one file in the AOI folder)\n",
    "for (dirpath, dirname, filename) in walk(aoi_folder):\n",
    "    aoi_filename = filename[0]\n",
    "    break\n",
    "    \n",
    "# Add file path\n",
    "aoi_path = aoi_folder + aoi_filename\n",
    "print(aoi_path) # Verify correct AOI file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "3c3fc3eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image CRS:\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "uint16 EPSG:32722 4\n",
      "-------------------\n",
      "AOI CRS: epsg:4326\n",
      "['epsg:4326', 'epsg:4326', 'epsg:4326', 'epsg:4326', 'epsg:4326', 'epsg:4326']\n"
     ]
    }
   ],
   "source": [
    "# Loading CRS of original images and AOI into corresponding lists\n",
    "img_list_crs = []\n",
    "aoi_list_crs = [0]*n_images\n",
    "aoi_info = fiona.open(aoi_path)\n",
    "    \n",
    "# Checking CRS of images and AOI\n",
    "print(\"Image CRS:\")\n",
    "for i in range(len(img_list)):\n",
    "    img_list_crs.append(img_list[i].meta['crs'])\n",
    "    aoi_list_crs[i] = aoi_info.crs['init']\n",
    "    print(img_list[i].meta['dtype'], img_list[i].meta['crs'], img_list[i].meta['count'])\n",
    "\n",
    "print('-------------------\\n' \"AOI CRS: {}\".format(aoi_list_crs[0]))\n",
    "print(aoi_list_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "97c93279",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create matching AOI CRS for each image CRS\n",
    "aoi_crs_transformed = []\n",
    "\n",
    "for i in range(len(img_list)):\n",
    "    # Use fiona to open the original AOI GeoJSON\n",
    "    with fiona.open(aoi_path) as mt:\n",
    "        aoi = [feature[\"geometry\"] for feature in mt]\n",
    "\n",
    "    # Transfrom AOI CRS to match image CRS\n",
    "    from fiona.transform import transform_geom\n",
    "    transformed_coords = transform_geom(str(aoi_list_crs[i]), str(img_list_crs[i]), aoi[0])\n",
    "    aoi = [transformed_coords]\n",
    "    aoi_crs_transformed.append([transformed_coords])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6924587",
   "metadata": {},
   "source": [
    "At this stage the AOI geometry is read and the coordinates have been transformed to match the images.\n",
    "Next, *rasterio.mask.mask* will create a mask over our images, using the AOI geometry as the mask line.\n",
    "\n",
    "Passing *crop=True* to the mask function will automatically crop the bits of our images that fall outside the mask boundary."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62399d18",
   "metadata": {},
   "source": [
    "#### Step 8. Export results\n",
    "Lastly, the result is saved as a final output GeoTIFF, scaled and clipped to selected AOI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ec86b650",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import rasterio's mask tool\n",
    "from rasterio.mask import mask\n",
    "\n",
    "# Loading clipped images into list\n",
    "img_list_clipped = [0]*n_images\n",
    "\n",
    "# Export reflectance images, scaled and clipped to AOI\n",
    "def scale_clip_export(images):\n",
    "    for i in range(len(images)):\n",
    "        \n",
    "        with rasterio.open(Results_path_scaled[i]) as image:\n",
    "            clipped, transform = mask(image, aoi, crop=True)\n",
    "            img_list_clipped[i] = clipped\n",
    "        \n",
    "        # Use metadata from our original images\n",
    "        meta = images[i].meta.copy()\n",
    "\n",
    "        # Update metadata with new, clipped image boundaries\n",
    "        meta.update({\"transform\": transform,\n",
    "            \"height\":clipped.shape[1],\n",
    "            \"width\":clipped.shape[2]})\n",
    "        \n",
    "        # Set the source metadata as kwargs we'll use to write the new data:\n",
    "        # update the 'dtype' value from uint16 to float64:\n",
    "        kwargs = meta\n",
    "        kwargs.update(dtype='uint16')\n",
    "        \n",
    "        # New image name for export\n",
    "        img_name_ref_scale_clip = (Results_folder + AnalyticMS_files[i]).replace('.tif','_Ref_clip.tif')\n",
    "\n",
    "        # Writes the new images into project folder\n",
    "        with rasterio.open(img_name_ref_scale_clip,'w', **kwargs) as dst:\n",
    "            dst.write(clipped)\n",
    "            \n",
    "        if i == (len(images)-1):\n",
    "                print(\"Success, {} scaled clipped images(dtype={}) have been exported to \"\\\n",
    "                \"your Results folder: {}\".format(n_images, kwargs['dtype'], Results_folder))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b11c0ec7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading clipped images into list\n",
    "img_list_clipped = [0]*n_images\n",
    "\n",
    "# Export reflectance images, unscaled and clipped to AOI\n",
    "def unscale_clip_export(images):\n",
    "    for i in range(len(images)):\n",
    "        \n",
    "        with rasterio.open(Results_path_unscaled[i]) as image:\n",
    "            clipped, transform = mask(image, aoi, crop=True)\n",
    "            img_list_clipped[i] = clipped\n",
    "        \n",
    "        # Use metadata from our original images\n",
    "        meta = images[i].meta.copy()\n",
    "\n",
    "        # Update metadata with new, clipped image boundaries\n",
    "        meta.update({\"transform\": transform,\n",
    "            \"height\":clipped.shape[1],\n",
    "            \"width\":clipped.shape[2]})\n",
    "        \n",
    "        # Set the source metadata as kwargs we'll use to write the new data:\n",
    "        # update the 'dtype' value from uint16 to float64:\n",
    "        kwargs = meta\n",
    "        kwargs.update(dtype='float64')\n",
    "        \n",
    "        # New image name for export\n",
    "        img_name_ref_scale_clip = (Results_folder + AnalyticMS_files[i]).\\\n",
    "        replace('.tif','_Ref_clip_unscaled.tif')\n",
    "\n",
    "        # Writes the new images into project folder\n",
    "        with rasterio.open(img_name_ref_scale_clip,'w', **kwargs) as dst:\n",
    "            dst.write(clipped)\n",
    "            \n",
    "        if i == (len(images)-1):\n",
    "                print(\"Success, {} unscaled clipped images(dtype={}) have been exported to \"\\\n",
    "                \"your Results folder: {}\".format(n_images, kwargs['dtype'], Results_folder))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ac26da3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Watch the results\n",
    "#show(img_list_clipped[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "755ac070",
   "metadata": {},
   "source": [
    "##### Export clipped Reflectance Images to your project folder\n",
    "- For the scaled clipped images, run the function *scale_clip_export()* with input(*img_list*)   \n",
    "\n",
    "\n",
    "- For the unscaled clipped images, run the function *unscale_clip_export()* with input(*img_list*)\n",
    "\n",
    "Note that unscaled(*dtype=float64*) takes 4 times as much storage space as scaled(*dtype=uint16*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "13f06ffb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success, 6 scaled clipped images(dtype=uint16) have been exported to your Results folder: planet_21-26/Results/\n"
     ]
    }
   ],
   "source": [
    "# Uncomment and run for wanted image type:\n",
    "# Export scaled images: - (default setting)\n",
    "scale_clip_export(img_list)\n",
    "\n",
    "# Export unscaled images:\n",
    "#unscale_clip_export(img_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adedab8f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "084b2d7c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}