Batch Processing v2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "20e37d9b",
   "metadata": {},
   "source": [
    "# Batch Processing Example\n",
    "In this example, we use the `micasense.imageset` class to load a set of directories of images into a list of `micasense.capture` objects, and we iterate over that list, saving out each image as an aligned stack of images as separate bands in a single tiff file each. Part of this process (via `imageutils.write_exif_to_stack`) injects that the GPS, capture datetime, camera model, etc into the processed images, allowing us to stitch those images using commercial software such as Pix4DMapper or Agisoft Metashape.\n",
    "\n",
    "Note: for this example to work, the images must have a valid RigRelatives tag. This requires RedEdge (3/M/MX) version of at least 3.4.0, or any version of RedEdge-P/Altum-PT/Altum/RedEdge-MX Dual. If your images don't meet that spec, you can also follow this support article to add the RigRelatives tag to your imagery: https://support.micasense.com/hc/en-us/articles/360006368574-Modifying-older-collections-for-Pix4Dfields-support"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c64ead96",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80646655",
   "metadata": {},
   "source": [
    "# Load Images into ImageSet\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2081b06b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import FloatProgress, Layout\n",
    "from IPython.display import display\n",
    "import micasense.imageset as imageset\n",
    "import micasense.capture as capture\n",
    "import os, glob\n",
    "import multiprocessing\n",
    "from pathlib import Path\n",
    "\n",
    "# set to True if you have an Altum-PT\n",
    "# or RedEdge-P and wish to output pan-sharpened stacks \n",
    "panSharpen = True \n",
    "\n",
    "# If creating a lot of stacks, it is more efficient to save the metadata\n",
    "# and then write all of the exif to the images after the stacks are created\n",
    "write_exif_to_individual_stacks = False\n",
    "\n",
    "panelNames = None\n",
    "useDLS = True\n",
    "\n",
    "# set your image path here. See more here: https://docs.python.org/3/library/pathlib.html\n",
    "imagePath = Path(\"./data/REDEDGE-MX\")\n",
    "\n",
    "# these will return lists of image paths as strings. Comment out of you aren't using panels. \n",
    "panelNames = list(imagePath.glob('IMG_0001_*.tif'))\n",
    "panelNames = [x.as_posix() for x in panelNames]\n",
    "\n",
    "if panelNames:\n",
    "    panelCap = capture.Capture.from_filelist(panelNames)\n",
    "\n",
    "# destinations on your computer to put the stacks\n",
    "# and RGB thumbnails\n",
    "outputPath = imagePath / '..' / 'stacks'\n",
    "thumbnailPath = outputPath / 'thumbnails'\n",
    "\n",
    "cam_model = panelCap.camera_model\n",
    "cam_serial = panelCap.camera_serial\n",
    "\n",
    "# determine if this sensor has a panchromatic band \n",
    "if cam_model == 'RedEdge-P' or cam_model == 'Altum-PT':\n",
    "    panchroCam = True\n",
    "else:\n",
    "    panchroCam = False\n",
    "    panSharpen = False \n",
    "    \n",
    "# if this is a multicamera system like the RedEdge-MX Dual,\n",
    "# we can combine the two serial numbers to help identify \n",
    "# this camera system later. \n",
    "if len(panelCap.camera_serials) > 1:\n",
    "    cam_serial = \"_\".join(panelCap.camera_serials)\n",
    "    print(\"Serial number:\",cam_serial)\n",
    "else:\n",
    "    cam_serial = panelCap.camera_serial\n",
    "    print(\"Serial number:\",cam_serial)\n",
    "    \n",
    "overwrite = False # can be set to set to False to continue interrupted processing\n",
    "generateThumbnails = True\n",
    "\n",
    "# Allow this code to align both radiance and reflectance images; but excluding\n",
    "# a definition for panelNames above, radiance images will be used\n",
    "# For panel images, efforts will be made to automatically extract the panel information\n",
    "# but if the panel/firmware is before Altum 1.3.5, RedEdge 5.1.7 the panel reflectance\n",
    "# will need to be set in the panel_reflectance_by_band variable.\n",
    "# Note: radiance images will not be used to properly create NDVI/NDRE images below.\n",
    "if panelNames is not None:\n",
    "    panelCap = capture.Capture.from_filelist(panelNames)\n",
    "else:\n",
    "    panelCap = None\n",
    "\n",
    "if panelCap is not None:\n",
    "    if panelCap.panel_albedo() is not None and not any(v is None for v in panelCap.panel_albedo()):\n",
    "        panel_reflectance_by_band = panelCap.panel_albedo()\n",
    "    else:\n",
    "        panel_reflectance_by_band = [0.49]*len(panelCap.eo_band_names()) #RedEdge band_index order\n",
    "    \n",
    "    panel_irradiance = panelCap.panel_irradiance(panel_reflectance_by_band)    \n",
    "    img_type = \"reflectance\"\n",
    "else:\n",
    "    if useDLS:\n",
    "        img_type='reflectance'\n",
    "    else:\n",
    "        img_type = \"radiance\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60736009",
   "metadata": {},
   "outputs": [],
   "source": [
    "## This progress widget is used for display of the long-running process\n",
    "f = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Loading\")\n",
    "display(f)\n",
    "def update_f(val):\n",
    "    if (val - f.value) > 0.005 or val == 1: #reduces cpu usage from updating the progressbar by 10x\n",
    "        f.value=val\n",
    "\n",
    "%time imgset = imageset.ImageSet.from_directory(imagePath, progress_callback=update_f)\n",
    "update_f(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c0a234d",
   "metadata": {},
   "source": [
    "# Capture map\n",
    "We can map out the capture GPS locations to ensure we are processing the right data. A GeoJSON of the captures will later be saved to the outputPath."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e9c437c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "from mapboxgl.viz import *\n",
    "from mapboxgl.utils import df_to_geojson, create_radius_stops, scale_between\n",
    "from mapboxgl.utils import create_color_stops\n",
    "import pandas as pd\n",
    "\n",
    "data, columns = imgset.as_nested_lists()\n",
    "df = pd.DataFrame.from_records(data, index='timestamp', columns=columns)\n",
    "\n",
    "#Insert your mapbox token here\n",
    "token = 'pk.eyJ1Ijoic3RlcGhlbm1hbmd1bTIiLCJhIjoiY2xmOXdnYzF1MDFqejNvdGE0YW13aTN5ZyJ9.AG_ckhUqTBjuGC2LuWCfQQ'\n",
    "color_property = 'dls-yaw'\n",
    "num_color_classes = 8\n",
    "\n",
    "min_val = df[color_property].min()\n",
    "max_val = df[color_property].max()\n",
    "\n",
    "import jenkspy\n",
    "geojson_data = df_to_geojson(df,columns[3:],lat='latitude',lon='longitude')\n",
    "breaks = jenkspy.jenks_breaks(df[color_property], nb_class=num_color_classes)\n",
    "color_stops = create_color_stops(breaks,colors='YlOrRd')\n",
    "\n",
    "viz = CircleViz(geojson_data, access_token=token, color_property=color_property,\n",
    "                color_stops=color_stops,\n",
    "                center=[df['longitude'].median(),df['latitude'].median()], \n",
    "                zoom=16, height='600px',\n",
    "                style='mapbox://styles/mapbox/satellite-streets-v9')\n",
    "viz.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e540a655",
   "metadata": {},
   "source": [
    "# Define which warp method to use\n",
    "For newer data sets with RigRelatives tags (images captured with RedEdge (3/M/MX) version 3.4.0 or greater with a valid calibration load, see https://support.micasense.com/hc/en-us/articles/360005428953-Updating-RedEdge-for-Pix4Dfields), we can use the RigRelatives for a simple alignment. To use this simple alignment, simply set `warp_matrices=None` \n",
    "\n",
    "For sets without those tags, or sets that require a RigRelatives optimization, we can go through the Alignment.ipynb notebook and get a set of warp_matrices that we can use here to align."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "635ad8d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array\n",
    "from numpy import float32\n",
    "from skimage.transform import ProjectiveTransform\n",
    "\n",
    "if panchroCam:\n",
    "    warp_matrices_filename = cam_serial + \"_warp_matrices_SIFT.npy\"\n",
    "else:\n",
    "    warp_matrices_filename = cam_serial + \"_warp_matrices_opencv.npy\"\n",
    "\n",
    "if Path('./' + warp_matrices_filename).is_file():\n",
    "    print(\"Found existing warp matrices for camera\", cam_serial)\n",
    "    load_warp_matrices = np.load(warp_matrices_filename, allow_pickle=True)\n",
    "    loaded_warp_matrices = []\n",
    "    for matrix in load_warp_matrices: \n",
    "        if panchroCam:\n",
    "            transform = ProjectiveTransform(matrix=matrix.astype('float64'))\n",
    "            loaded_warp_matrices.append(transform)\n",
    "        else:\n",
    "            loaded_warp_matrices.append(matrix.astype('float32'))\n",
    "\n",
    "    if panchroCam:\n",
    "        warp_matrices_SIFT = loaded_warp_matrices\n",
    "    else:\n",
    "        warp_matrices = loaded_warp_matrices\n",
    "    print(\"Loaded warp matrices from\",Path('./' + warp_matrices_filename).resolve())\n",
    "else:\n",
    "    print(\"No warp matrices found at expected location:\",warp_matrices_filename)\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aae7f09",
   "metadata": {},
   "source": [
    "## Align images and save each capture to a layered TIFF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc195da",
   "metadata": {},
   "outputs": [],
   "source": [
    "import exiftool\n",
    "import datetime\n",
    "import micasense.imageutils as imageutils\n",
    "exif_list = []\n",
    "## This progress widget is used for display of the long-running process\n",
    "f2 = FloatProgress(min=0, max=1, layout=Layout(width='100%'), description=\"Saving\")\n",
    "display(f2)\n",
    "def update_f2(val):\n",
    "    f2.value=val\n",
    "\n",
    "if not os.path.exists(outputPath):\n",
    "    os.makedirs(outputPath)\n",
    "if generateThumbnails and not os.path.exists(thumbnailPath):\n",
    "    os.makedirs(thumbnailPath)\n",
    "\n",
    "# Save out geojson data so we can open the image capture locations in our GIS\n",
    "with open(os.path.join(outputPath,'imageSet.json'),'w') as f:\n",
    "    f.write(str(geojson_data))\n",
    "    \n",
    "try:\n",
    "    irradiance = panel_irradiance+[0]\n",
    "except NameError:\n",
    "    irradiance = None\n",
    "\n",
    "start = datetime.datetime.now()\n",
    "for i,capture in enumerate(imgset.captures):\n",
    "    outputFilename = str(i).zfill(4) + \"_\" + capture.uuid+'.tif'\n",
    "    thumbnailFilename = str(i).zfill(4) + \"_\" + capture.uuid+'.jpg'\n",
    "    fullOutputPath = os.path.join(outputPath, outputFilename)\n",
    "    fullThumbnailPath= os.path.join(thumbnailPath, thumbnailFilename)\n",
    "    if (not os.path.exists(fullOutputPath)) or overwrite:\n",
    "        if(len(capture.images) == len(imgset.captures[0].images)):\n",
    "            if panchroCam:\n",
    "                capture.radiometric_pan_sharpened_aligned_capture(warp_matrices=warp_matrices_SIFT,irradiance_list=capture.dls_irradiace(), img_type=img_type, write_exif=write_exif_to_individual_stacks)\n",
    "            else:\n",
    "                capture.create_aligned_capture(irradiance_list=irradiance, warp_matrices=warp_matrices)\n",
    "            exif_list.append(imageutils.prepare_exif_for_stacks(capture,fullOutputPath))\n",
    "            capture.save_capture_as_stack(fullOutputPath, pansharpen=panSharpen,sort_by_wavelength=True, write_exif=write_exif_to_individual_stacks)\n",
    "            if generateThumbnails:\n",
    "                capture.save_capture_as_rgb(fullThumbnailPath)\n",
    "    capture.clear_image_data()\n",
    "    update_f2(float(i)/float(len(imgset.captures)))\n",
    "update_f2(1.0)\n",
    "end = datetime.datetime.now()\n",
    "\n",
    "print(\"Saving time: {}\".format(end-start))\n",
    "print(\"Alignment+Saving rate: {:.2f} images per second\".format(float(len(imgset.captures))/float((end-start).total_seconds())))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4b88c80",
   "metadata": {},
   "source": [
    "# Write EXIF data to stacks\n",
    "As mentioned above, it is more time intensive to write the exif data to each image as it is created. Here, we write the exif data after all of the TIFF files have been created. This should take a few seconds per stack."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac88b6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if write_exif_to_individual_stacks == False:\n",
    "    start = datetime.datetime.now()\n",
    "    for exif in exif_list:\n",
    "        imageutils.write_exif_to_stack(existing_exif_list=exif)\n",
    "    end = datetime.datetime.now()\n",
    "    print(\"Saving time: {}\".format(end-start))\n",
    "    print(\"Alignment+Saving rate: {:.2f} images per second\".format(float(len(exif_list))/float((end-start).total_seconds())))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea464070",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}