PS_preprocess_reflectance_mosaic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "52ce4a11",
   "metadata": {},
   "source": [
    "### Planetscope preprocessing - Reflectance Mosaic\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46a9410c",
   "metadata": {},
   "source": [
    "This script:\n",
    "- Reads multiple Planetscope images/scenes\n",
    "- Convert TOA Radiance to TOA Reflectance\n",
    "- Creates a composite image (or mosaic) \n",
    "- Clips mosaic 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. June 2019 over a selected study area in the Brazilian Amazon.\n",
    "\n",
    "Available export options:\n",
    "\n",
    "- Reflectance images scaled - *(default setting)*\n",
    "- Reflectance mosaic scaled - *(default setting)*\n",
    "- Reflectance mosaic scaled clipped to AOI - *(default setting)*   \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "- Reflectance images unscaled\n",
    "- Reflectance mosaic unscaled\n",
    "- Reflectance mosaic unscaled clipped to AOI\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a6407d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import rasterio\n",
    "import fiona"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "066a9319",
   "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 reflectance images, mosaic and clipped mosaic.\n",
    "\n",
    "The unscaled images, mosaic and unscaled clipped mosaic are available by calling the *unscaled-functions* in Step 6, 7 and 9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1bd0f947",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_images = 7\n",
    "Project_folder = 'planet_21/'\n",
    "AnalyticMS_folder = 'planet_21/AnalyticMS_files/'\n",
    "Metadata_folder = 'planet_21/Metadata_files/'\n",
    "aoi_folder = 'planet_21/AOI/'\n",
    "Results_folder = 'planet_21/Results/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cea858a7",
   "metadata": {},
   "source": [
    "#### Step 2. Read and Inspect Images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbe3632c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "planet_21/AnalyticMS_files/20190721_133241_1034_3B_AnalyticMS.tif\n",
      "planet_21/AnalyticMS_files/20190721_133242_1034_3B_AnalyticMS.tif\n",
      "planet_21/AnalyticMS_files/20190721_133243_1034_3B_AnalyticMS.tif\n",
      "planet_21/AnalyticMS_files/20190721_133316_1044_3B_AnalyticMS.tif\n",
      "planet_21/AnalyticMS_files/20190721_133317_1044_3B_AnalyticMS.tif\n",
      "planet_21/AnalyticMS_files/20190721_133318_1044_3B_AnalyticMS.tif\n",
      "planet_21/AnalyticMS_files/20190721_133319_1044_3B_AnalyticMS.tif\n"
     ]
    }
   ],
   "source": [
    "# Use 7 scenes from 21. June 2019 for coversion to TOA Reflectance mosaic\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": "64c9f941",
   "metadata": {},
   "source": [
    ""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "505a2e22",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "planet_21/Metadata_files/20190721_133241_1034_3B_AnalyticMS_metadata.xml\n",
      "planet_21/Metadata_files/20190721_133242_1034_3B_AnalyticMS_metadata.xml\n",
      "planet_21/Metadata_files/20190721_133243_1034_3B_AnalyticMS_metadata.xml\n",
      "planet_21/Metadata_files/20190721_133316_1044_3B_AnalyticMS_metadata.xml\n",
      "planet_21/Metadata_files/20190721_133317_1044_3B_AnalyticMS_metadata.xml\n",
      "planet_21/Metadata_files/20190721_133318_1044_3B_AnalyticMS_metadata.xml\n",
      "planet_21/Metadata_files/20190721_133319_1044_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: 7 \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": "markdown",
   "id": "f747cc08",
   "metadata": {},
   "source": [
    "At this point we can use rasterio to inspect the metadata of these three images. Specifically, in order to create a composite from these images, we want to verify that all three images have the same data type, the same coordinate reference systems and the same band count"
   ]
  },
  {
   "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",
      "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.\n",
    "\n",
    "Now that we've verified all three satellite images have the same critical metadata, we can safely use rasterio.merge to stitch them together."
   ]
  },
  {
   "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.22787121429e-05, 2: 2.3597734536e-05, 3: 2.62801712705e-05, 4: 3.9553550288e-05}\n",
      "{1: 2.22957954184e-05, 2: 2.36158292355e-05, 3: 2.63003228576e-05, 4: 3.95838798777e-05}\n",
      "{1: 2.23129263982e-05, 2: 2.36339744636e-05, 3: 2.6320530717e-05, 4: 3.96142941613e-05}\n",
      "{1: 2.22858931046e-05, 2: 2.36138635285e-05, 3: 2.63757407995e-05, 4: 3.9632589334e-05}\n",
      "{1: 2.23029452396e-05, 2: 2.3631931765e-05, 3: 2.63959222968e-05, 4: 3.96629143589e-05}\n",
      "{1: 2.23200437529e-05, 2: 2.36500491436e-05, 3: 2.64161586838e-05, 4: 3.96933218618e-05}\n",
      "{1: 2.23372065248e-05, 2: 2.36682346096e-05, 3: 2.64364711218e-05, 4: 3.97238436401e-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-5, 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",
    "\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 10839\n",
      "Red band reflectance goes from: 0.0 to 0.28485077640094947\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": "42f407bc",
   "metadata": {},
   "source": [
    "#### Step 6. Write the Reflectance Images\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": "ca7957fa",
   "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": "0dbcb3f4",
   "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": "e009d03f",
   "metadata": {},
   "source": [
    "Export of unscaled GeoTIFF, dtype=float64 with Reflectance values between 0.0 - 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5704a3a7",
   "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 from uint16 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": "e4775e53",
   "metadata": {},
   "source": [
    "##### Export Reflectance Images to your project 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 mosaicing and clipping to AOI, proceed further"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2c0726db",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With scale factor: 10000\n",
      "Before scaling: Red band reflectance from: 0.0 to 0.28485077640094947\n",
      "After scaling: Red band reflectance from: 0.0 to 2848.5077640094946\n",
      "Success, 7 scaled images(dtype=uint16) have been exported to your Results folder: planet_21/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": "80a79de7",
   "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": "f82ed8c7",
   "metadata": {},
   "source": [
    "Using rasterio.plot (a matplotlib interface) to preview the results of our mosaic."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1358f3b",
   "metadata": {},
   "source": [
    "#### Step 7. Creating the Mosaic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5a43d04c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge returns the mosaic & coordinate transformation information\n",
    "from rasterio.merge import merge\n",
    "(mosaic_scaled, transform) = merge(img_list_scaled)\n",
    "(mosaic_unscaled, transform) = merge(img_list_unscaled)"
   ]
  },
  {
   "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": 29,
   "id": "311ed615",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preview results if desired\n",
    "from rasterio.plot import show\n",
    "#show(mosaic_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4688b7a4",
   "metadata": {},
   "source": [
    "Before writing the mosaic out to a new GeoTIFF file, the metadata is taken from the first image/scene to represent the metadata of all input images.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "808c2b36",
   "metadata": {},
   "source": [
    "Export of scaled unclipped mosaic GeoTIFF (multiplied with scale factor of 10 000 and stored as *dtype=uint16*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cd5b95e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to export the unclipped scaled mosaic\n",
    "def mosaic_scaled_export(images):\n",
    "    # Get the metadata of original GeoTIFF:\n",
    "        meta = images[0].meta\n",
    "        \n",
    "        # Update the original metadata to reflect the specifics of our new mosaic\n",
    "        meta.update({\"transform\": transform,\n",
    "                     \"height\":mosaic_scaled.shape[1],\n",
    "                     \"width\":mosaic_scaled.shape[2]})\n",
    "        \n",
    "        # New name for exported scaled mosaic\n",
    "        img_name_mosaic_scaled = (Results_folder + 'Ref_Mosaic.tif')\n",
    "        \n",
    "        # Write band calculations to a new scaled GeoTIFF file\n",
    "        with rasterio.open(img_name_mosaic_scaled, 'w', **meta) as dst:\n",
    "                dst.write(mosaic_scaled)\n",
    "\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, a scaled mosaic(dtype={}) composed of {} images have been \"\\\n",
    "                  \"exported to your Results folder: {}\".format(meta['dtype'], n_images, Results_folder ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c3a407f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to export the unclipped scaled mosaic\n",
    "def mosaic_unscaled_export(images):\n",
    "    # Get the metadata of original GeoTIFF:\n",
    "        meta = images[0].meta\n",
    "        \n",
    "        # Update the original metadata to reflect the specifics of our new mosaic\n",
    "        meta.update({\"transform\": transform,\n",
    "                     \"height\":mosaic_unscaled.shape[1],\n",
    "                     \"width\":mosaic_unscaled.shape[2]})\n",
    "        \n",
    "        # update the 'dtype' value to float64:\n",
    "        kwargs = meta\n",
    "        meta.update(dtype='float64')\n",
    "        \n",
    "        # New name for exported mosaic unscaled\n",
    "        img_name_mosaic_unscaled = (Results_folder + 'Ref_Mosaic_unscaled.tif')\n",
    "        \n",
    "        # Write band calculations to a new unscaled GeoTIFF file\n",
    "        with rasterio.open(img_name_mosaic_unscaled, 'w', **meta) as dst:\n",
    "                dst.write(mosaic_unscaled)\n",
    "\n",
    "        if i == (len(images)-1):\n",
    "            print(\"Success, a unscaled mosaic(dtype={}) composed of {} images have been \"\\\n",
    "                  \"exported to your Results folder: {}\".format(meta['dtype'], n_images, Results_folder ))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8f5b78",
   "metadata": {},
   "source": [
    "##### Export Reflectance Mosaic to your Results folder\n",
    "- For the scaled mosaic, run the function *mosaic_scaled_export()* with input(*img_list*)   \n",
    "\n",
    "\n",
    "- For the unscaled images, run the function *mosaic_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 clipped reflectance mosaic to AOI, proceed further"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "8ca548f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success, a scaled mosaic(dtype=uint16) composed of 7 images have been exported to your Results folder: planet_21/Results/\n"
     ]
    }
   ],
   "source": [
    "# Uncomment and run for wanted mosaic type:\n",
    "# Export scaled mosaic:\n",
    "mosaic_scaled_export(img_list)\n",
    "\n",
    "# Export unscaled mosaic:\n",
    "#mosaic_unscaled_export(img_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dfc6500",
   "metadata": {},
   "source": [
    "#### Step 8. Clip the Mosaic to AOI Boundaries\n",
    "A mask is used to clip the mosaic to the extents of our area of interest (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 mosaic, 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": 21,
   "id": "b1bea761",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "planet_21/AOI/LCC_AOI.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": 22,
   "id": "b3d2696c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mosaic CRS: EPSG:32722\n",
      "-------------------\n",
      "AOI CRS: epsg:4326\n"
     ]
    }
   ],
   "source": [
    "# Checking CRS of images and AOI\n",
    "aoi_info = fiona.open(aoi_path)\n",
    "mosaic_crs = img_list[0].meta['crs']\n",
    "aoi_crs = aoi_info.crs['init']\n",
    "    \n",
    "print(\"Mosaic CRS:\", mosaic_crs)\n",
    "print('-------------------\\n' \"AOI CRS: {}\".format(aoi_crs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4b9e7286",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 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 mosaic CRS\n",
    "from fiona.transform import transform_geom\n",
    "transformed_coords = transform_geom(str(aoi_crs), str(mosaic_crs), aoi[0])\n",
    "aoi = [transformed_coords]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b2026acf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import rasterio's mask tool\n",
    "from rasterio.mask import mask\n",
    "\n",
    "# apply mask with crop=True to cut to boundary\n",
    "with rasterio.open('planet_21/Results/Ref_Mosaic.tif') as mosaic:\n",
    "    clipped_scaled, transform = mask(mosaic, aoi, crop=True)\n",
    "    \n",
    "with rasterio.open('planet_21/Results/Ref_Mosaic_unscaled.tif') as mosaic:\n",
    "    clipped_unscaled, transform = mask(mosaic, aoi, crop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "150b6f43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Watch results\n",
    "#show(clipped_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cef51fb2",
   "metadata": {},
   "source": [
    "#### Step 9. Export results\n",
    "Lastly, the result is saved as a final output GeoTIFF, scaled and clipped to selected AOI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9c875161",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the output to a final GeoTIFF\n",
    "def mosaic_scale_clip_export(images):\n",
    "    # use the metadata from our original mosaic\n",
    "    meta = images[0].meta.copy()\n",
    "\n",
    "    # update metadata with new, clipped mosaic's boundaries\n",
    "    meta.update({\"transform\": transform,\n",
    "        \"height\":clipped_scaled.shape[1],\n",
    "        \"width\":clipped_scaled.shape[2]})\n",
    "\n",
    "    # write the clipped scaled mosaic to a GeoTIFF\n",
    "    with rasterio.open(Results_folder + \"Ref_Mosaic_Clip.tif\", 'w', **meta) as dst:\n",
    "        dst.write(clipped_scaled)\n",
    "    \n",
    "    print(\"Success, a scaled clipped mosaic(dtype={}) composed of {} images have been \"\\\n",
    "    \"exported to your Results folder: {}\".format(meta['dtype'], n_images, Results_folder ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "c4b3ecfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def mosaic_unscale_clip_export(images):\n",
    "    # use the metadata from our original mosaic\n",
    "    meta = img_list_unscaled[0].meta.copy()\n",
    "\n",
    "    # update metadata with new, clipped mosaic's boundaries\n",
    "    meta.update({\"transform\": transform,\n",
    "        \"height\":clipped_scaled.shape[1],\n",
    "        \"width\":clipped_scaled.shape[2]})\n",
    "\n",
    "    # update the 'dtype' value to float64:\n",
    "    kwargs = meta\n",
    "    kwargs.update(dtype='float64')\n",
    "\n",
    "    # write the clipped unscaled mosaic to a GeoTIFF\n",
    "    with rasterio.open(Results_folder + \"Ref_Mosaic_Clip_unscaled.tif\", 'w', **kwargs) as dst:\n",
    "        dst.write(clipped_unscaled)\n",
    "    \n",
    "    print(\"Success, a unscaled clipped mosaic(dtype={}) composed of {} images have been \"\\\n",
    "    \"exported to your Results folder: {}\".format(kwargs['dtype'], n_images, Results_folder ))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "781c981d",
   "metadata": {},
   "source": [
    "##### Export clipped Reflectance Mosaics to your Results folder\n",
    "- For the scaled mosaics, run the function *mosaic_scale_clip_export()* with input(*img_list*)   \n",
    "\n",
    "\n",
    "- For the unscaled mosaics, run the function *mosaic_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": 27,
   "id": "fe51c10f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success, a scaled clipped mosaic(dtype=uint16) composed of 7 images have been exported to your Results folder: planet_21/Results/\n"
     ]
    }
   ],
   "source": [
    "# Uncomment and run for wanted mosaic type:\n",
    "# Export scaled clipped mosaic:\n",
    "mosaic_scale_clip_export(img_list)\n",
    "\n",
    "# Export unscaled clipped mosaic:\n",
    "#mosaic_unscale_clip_export(img_list)"
   ]
  }
 ],
 "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
}