introduction.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyOgVDQ2GSRiRB8pBVoDflDN",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/wurDevTim/Workshop_P4P/blob/main/introduction.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Welkom to the P4P workshop!\n",
        "This notebook gives a few examples on how to analyze cropreporter data with plantcv.\n",
        "\n",
        "There are also various tutorials on the plantcv website, unfortunately not all are updated to plantcv v4.0: https://plantcv.readthedocs.io/en/stable/tutorials/\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "-7ASH8IhV327"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Preperation\n",
        "For this training we prepared a small dataset you can download from:\n",
        "https://drive.google.com/file/d/1i-EIXZ-hw4uwKRV40rPAql9Lb-eIDm1j/view?usp=sharing\n",
        "\n",
        "Unzip the folder and upload the content to your own google drive."
      ],
      "metadata": {
        "id": "KXnETprgbxEp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Plantcv only suports python version 3.8, 3.9 and 3.10.\n",
        "# Colab only supports one python version, the current version is 3.10.12\n",
        "import sys\n",
        "print(f\"Running python version: {sys.version}\")"
      ],
      "metadata": {
        "id": "5_htLaMHYSuz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Install the requirements"
      ],
      "metadata": {
        "id": "JVogt2aWgjlj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# plantcv requires a newer version of altair than installed by default\n",
        "!pip install altair==5.0.1\n",
        "!pip install plantcv"
      ],
      "metadata": {
        "id": "77-MeqVjV9Iz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Mount google drive\n",
        "from google.colab import drive\n",
        "from os import path\n",
        "\n",
        "datafolder = \"/content/drive/My Drive/P4P_workshop_data\"\n",
        "# Check if the data folder is mounted correctly\n",
        "if not path.exists(datafolder):\n",
        "  drive.mount('/content/drive')\n",
        "\n",
        "!ls \"$datafolder\""
      ],
      "metadata": {
        "id": "oUpDB9x64mqC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Importing Libraries\n",
        "from google.colab.patches import cv2_imshow\n",
        "# Import packages\n",
        "from plantcv import plantcv as pcv\n",
        "from skimage.util import img_as_ubyte\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import cv2"
      ],
      "metadata": {
        "id": "d3vZXSIWwkRy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "print(f\"Currently installed plantcv version is v{pcv.__version__}\")"
      ],
      "metadata": {
        "id": "AGv-gwJhoM0S"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Plantcv\n",
        "Plantcv implements quite some functions, in the notebook examples are gives of how to load the data and:\n",
        "- Use the Chlorophyll image to create a mask\n",
        "- Compute the spectral indices\n",
        "- Compute the FV/FM\n",
        "- Analyse an entire folder automatically\n",
        "With two bonusses at the end\n",
        "\n",
        "For more information see:\n",
        "https://plantcv.readthedocs.io/en/stable/photosynthesis_read_cropreporter/"
      ],
      "metadata": {
        "id": "ZM4Vw_PFkESF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Configure plantcv to show debug output\n",
        "pcv.params.debug = \"plot\"\n",
        "pcv.params.debug_outdir = \"temp_output\"\n",
        "outdir = 'temp_output'"
      ],
      "metadata": {
        "id": "6tTUdZOg3Qry"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Plantcv photosynthesis object\n",
        "To start using plantcv for analysis of cropreporter data the following function can be called with the `.inf` file:\n",
        "\n",
        "```\n",
        "pcv.photosynthesis.read_cropreporter\n",
        "```\n",
        "The function automatically checks the folder for other cropreporter files like the SPC, PML, PMD, CHL and CLR file.\n",
        "\n",
        "The function returns a PSII_data instance, which contains a subclass per data object:\n",
        "*   ojip_light\n",
        "*   ojip_dark\n",
        "*   pam_light\n",
        "*   pam_dark\n",
        "*   chlorophyll\n",
        "*   spectral\n",
        "\n",
        "The RGB data is collected with the same monochrome camera as the spectral data, the only difference is the filter that is used: a Red, green or Blue filter.\n",
        "It can be accessed with: `ps.spectral.pseudo_rgb`.\n",
        "If a parameter was not measured the class is set to `None`."
      ],
      "metadata": {
        "id": "4MQP0SbM9bm0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Load the first dataset into plantcv\n",
        "\n",
        "# Note this dataset does not contain photosynthesis (PAM) data\n",
        "#filename = datafolder + \"/20231018_AO1/HDR_96_NPEC51.20230906.AO1.MoneyMaker.Salt.88_1476.INF\"\n",
        "\n",
        "# Note: This dataset does not contain spectral(SPC) data\n",
        "#filename = datafolder + \"/Enza Mandy Boon/ChiVMV/HDR_E0044P0206N0001.INF\"\n",
        "\n",
        "# Note: This dataset does not contain color or spectral (CLR or SPC) data\n",
        "#filename = datafolder + \"/Marie/NPEC51.20231124.AA12.CYCTF23_0002.CYC_23_Y_BT1E.6/CropReporter/181215797/data/HDR_100_NPEC51.20231124.AA12.CYCTF23_0002.CYC_23_Y_BT1E.6_1515.INF\"\n",
        "\n",
        "# Note: This dataset does not contain spectral(SPC) data\n",
        "#filename = datafolder + \"/zheng/resiatant materials_1/0hpi/HDR_Experiment20230802shangwu_B1_P27_U1_DIRXpY0_T_G_2023-08-09_23-01-23.INF\"\n",
        "\n",
        "# Contains most data types, therefore this is ideal to start with\n",
        "filename = datafolder + \"/Lucia/NPEC52.20230605.BD4.CE3027.Control.1/CropReporter/124326957/data/HDR_90_NPEC52.20230605.BD4.CE3027.Control.1_1368.INF\"\n",
        "\n",
        "# Load the data\n",
        "ps = pcv.photosynthesis.read_cropreporter(filename=filename)"
      ],
      "metadata": {
        "id": "_JVE-QAjV9g2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Show the chlorophyll image\n",
        "chlorophyll = ps.chlorophyll.sel(frame_label='Chl').to_numpy()\n",
        "plt.imshow(chlorophyll)"
      ],
      "metadata": {
        "id": "-PAfC74apoXS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Note the big difference between the plant and background.\n",
        "# this makes the chlorophyll ideall to create masks\n",
        "mask = pcv.threshold.otsu(gray_img=img_as_ubyte(ps.chlorophyll.sel(frame_label=\"Chl\").data), object_type=\"light\")\n",
        "#mask = pcv.fill_holes(bin_img=mask)"
      ],
      "metadata": {
        "id": "stoTtrq63EIf"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Compute average chlorophyll\n",
        "chlorophyll = ps.chlorophyll.sel(frame_label='Chl').to_numpy()\n",
        "masked_chlorophyll = cv2.bitwise_and(chlorophyll, chlorophyll, mask=mask)\n",
        "mean_chlorophyll = np.nanmean(chlorophyll[mask > 0])\n",
        "print(f\"average chlorophyll: {round(mean_chlorophyll, 2)}\")"
      ],
      "metadata": {
        "id": "0QLVEYzoAm6V"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Show the RGB image\n",
        "img = ps.spectral.pseudo_rgb\n",
        "img = cv2.bitwise_and(img, img, mask=mask)\n",
        "# first parameter is the contrast scaling, the second the brightness, which is increased slightly for visualization\n",
        "# For a concrete example see: https://www.tutorialspoint.com/how-to-change-the-contrast-and-brightness-of-an-image-using-opencv-in-python\n",
        "img = cv2.convertScaleAbs(img, 1.0, 10)\n",
        "\n",
        "f, axarr = plt.subplots(1,2, figsize=(15, 15))\n",
        "axarr[0].imshow(cv2.convertScaleAbs(ps.spectral.pseudo_rgb, 1.0, 10))\n",
        "axarr[1].imshow(img)\n",
        "plt.show()\n",
        "# Left the original image, Right the masked image"
      ],
      "metadata": {
        "id": "3scGbgA_-W_k"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "It is challenging to find a solution that scales to all plant species.\n",
        "Here will show a few corrections which can be applied.\n",
        "\n",
        "### Exercise\n",
        "Try loading a different file by using one of the different filenames and play with the masking. Would you always use the chlorophyll image? Mask on the Photosynthesis data or use something else?\n",
        "\n",
        "Below two examples are shown of how to improve the mask"
      ],
      "metadata": {
        "id": "0Pt_W0zfpuNg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Mask below chlorophyll value\n",
        "threshold = 4000\n",
        "chlorophyll = ps.chlorophyll.sel(frame_label='Chl').to_numpy()\n",
        "threshold_mask = np.zeros(chlorophyll.shape, dtype=np.uint8)\n",
        "idx = np.where(chlorophyll > threshold)\n",
        "threshold_mask[idx] = 255\n",
        "plt.imshow(threshold_mask, cmap='gray')"
      ],
      "metadata": {
        "id": "G5t8qOUmtLyA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Fill small objects to remove noise and get a complete plant\n",
        "# Inputs:\n",
        "#   bin_img         - Binary image data\n",
        "#   size            - Minimum object area size in pixels (integer), smaller objects get filled in.\n",
        "filled_mask = pcv.fill(bin_img=mask, size=200)"
      ],
      "metadata": {
        "id": "OId0au1cqPBb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Erode a small border of pixels from the mask\n",
        "# Inputs:\n",
        "#   gray_img - Grayscale (usually binary) image data\n",
        "#   ksize - The size used to build a ksize x ksize\n",
        "#            matrix using np.ones. Must be greater than 1 to have an effect\n",
        "#   i - An integer for the number of iterations\n",
        "eroded_mask = pcv.erode(gray_img=mask, ksize=3, i=1)"
      ],
      "metadata": {
        "id": "ZK9lb8OfqhD5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "unmasked_rgb = cv2.convertScaleAbs(ps.spectral.pseudo_rgb, 1.0, 10)\n",
        "f, axarr = plt.subplots(1,2, figsize=(15, 15))\n",
        "axarr[0].imshow(cv2.bitwise_and(unmasked_rgb, unmasked_rgb, mask=filled_mask))\n",
        "axarr[1].imshow(cv2.bitwise_and(unmasked_rgb, unmasked_rgb, mask=eroded_mask))\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "I5LeAtBEqqhh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# One way to compare masks on similairy is to calculate the intersection over union.\n",
        "# It returns values between 0-1, where 1 is identical and zero shows no overlap.\n",
        "def calculate_iou(maska, maskb):\n",
        "  maska = maska.flatten()\n",
        "  maskb = maskb.flatten()\n",
        "  mask = np.zeros_like(maska)\n",
        "  mask[maska > 0] = 1\n",
        "  mask_two = np.zeros_like(maskb)\n",
        "  mask_two[maskb > 0] = 1\n",
        "\n",
        "  overlap = mask * mask_two  # Logical AND\n",
        "  union = (mask + mask_two)>0  # Logical OR\n",
        "  iou = overlap.sum() / float(union.sum())\n",
        "  return iou\n",
        "\n",
        "print(f\"Similairy of the original mask and filled mask: {round(calculate_iou(mask, filled_mask), 3)}, similairy of the erorded and original mask: {round(calculate_iou(mask, eroded_mask), 3)}\")"
      ],
      "metadata": {
        "id": "fKOKMPmC403z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Spectral indices\n",
        "Different pigments and chlorophyll content can be measured at different wavelengths, by combining these measurements the relation can be analyzed. As example, NDVI visualizes changes in the chlorophyl absorption of the plant.\n",
        "A lower value indicates somethings is wrong with the plant, but it does not say anything about the cause.\n",
        "It is also important to note that:\n",
        "- Each camera uses slightly different settings/configurations which result in different NDVI values, making it challenging to compare or exactly reproduce experiments.\n",
        "- Plants grow towards the camera, as they come closer there is less dispersion of light which is no problem when you fly a drone 30 meters above a potato field but introduces bias in measurement taken up-close, as done by the cropreporter.\n",
        "\n",
        "There is an extensive list of spectral indexes available in plantcv, however not all can be computer with the wavelengths measured by the cropreporter:\n",
        "https://github.com/danforthcenter/plantcv/blob/main/plantcv/plantcv/spectral_index/spectral_index.py\n",
        "\n",
        "The indexes we show here are chosen based on this study:\n",
        "[Application of Phenotyping Methods in Detection of Drought and Salinity Stress in Basil (Ocimum basilicum L.) (Boris et al. 2021)](https://www.frontiersin.org/articles/10.3389/fpls.2021.629441/full)\n",
        "- NDVI\n",
        "- ARI\n",
        "- MARI\n",
        "- PSRI\n",
        "- Egreen"
      ],
      "metadata": {
        "id": "Rbi1Q01XckC7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Wavelengths\n",
        "A few weeks ago, we discovered a small bug in the plantcv code, which is described here [cropreporter wavelengths](https://github.com/danforthcenter/plantcv/issues/1373)\n",
        "This is fixed in the main branch, as there is no new release, we are currently using a plantcv version which contains this bug.\n",
        "\n",
        "Currently the wavelengths are: {460, 500, 550, 670, 700, 800}\n",
        "this be checked by running: `ps.spectral.wavelength_dict`\n",
        "\n",
        "The wavelengths of the NPEC cropreporter are in reality: {475, 540, 550, 640, 710, 770}\n",
        "Unfortunately, phenovation does not store this in their metadata so we cannot guaranty that this is the same for all systems sold worldwide.\n"
      ],
      "metadata": {
        "id": "2UiP8aFq0Q8g"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# In the fix below we asume the dataset contains both a CLR and SPC file.\n",
        "# Show the wavelengths\n",
        "ps.spectral.wavelength_dict"
      ],
      "metadata": {
        "id": "_XGVV4540Pq0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Overwrite the dictionary with the correct wavelengths\n",
        "ps.spectral.wavelength_dict = {475: 0.0, 540: 1.0, 550: 2.0, 640: 3.0, 710: 4.0, 770: 5.0}"
      ],
      "metadata": {
        "id": "sYPexUFt2Nlz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Normalized Difference Vegetation Index: (R800 - R670) / (R800 + R670)\n",
        "# The theoretical range for NDVI is [-1.0, 1.0]\n",
        "ndvi_img = pcv.spectral_index.ndvi(hsi=ps.spectral, distance=30).array_data\n",
        "masked_ndvi_img = cv2.bitwise_and(ndvi_img, ndvi_img, mask=mask)\n",
        "mean_ndvi = np.nanmean(ndvi_img[mask > 0])\n",
        "print(f\"ndvi: {mean_ndvi}\")\n",
        "plt.imshow(masked_ndvi_img)"
      ],
      "metadata": {
        "id": "_zeOlsb3dHBg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Anthocyanin Reflectance Index: (1 / R550) - (1 / R700)\n",
        "# The theoretical range for ARI is (-Inf, Inf)\n",
        "ari_img = pcv.spectral_index.ari(hsi=ps.spectral).array_data\n",
        "masked_ari_img = cv2.bitwise_and(ari_img, ari_img, mask=mask)\n",
        "plt.imshow(masked_ari_img)\n",
        "mean_ari = np.nanmean(ari_img[mask > 0])\n",
        "print(f\"ari: {mean_ari}\")"
      ],
      "metadata": {
        "id": "ZSuKwbFXdZCQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Modified Anthocyanin Reflectance Index: ((1 / R550) - (1 / R700)) * R800\n",
        "# The theoretical range for Modified ARI is (-Inf, Inf)\n",
        "mari_img = pcv.spectral_index.mari(hsi=ps.spectral, distance=30).array_data\n",
        "masked_mari_img = cv2.bitwise_and(mari_img, mari_img, mask=mask)\n",
        "plt.imshow(masked_mari_img)\n",
        "mean_mari = np.nanmean(mari_img[mask > 0])\n",
        "print(f\"modified ari: {mean_mari}\")"
      ],
      "metadata": {
        "id": "eCXqLqPEeDhg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plant Senescence Reflectance Index: (R678 - R500) / R750\n",
        "# The theoretical range for PSRI is (-Inf, Inf)\n",
        "# An increase in PSRI sugest an increase in canopy stress.\n",
        "psri_img = pcv.spectral_index.psri(hsi=ps.spectral).array_data\n",
        "masked_psri_img = cv2.bitwise_and(psri_img, psri_img, mask=mask)\n",
        "plt.imshow(masked_psri_img)\n",
        "mean_psri = np.nanmean(psri_img[mask > 0])\n",
        "print(f\"psri: {mean_psri}\")"
      ],
      "metadata": {
        "id": "1dx8SGoAdwR6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Excess Green Index: 2g - r - b\n",
        "# The theoretical range for EGI is (-1, 2)\n",
        "# Question why do you see warnings?\n",
        "egreen_img = pcv.spectral_index.egi(ps.spectral.pseudo_rgb).array_data\n",
        "masked_egreen_img = cv2.bitwise_and(egreen_img, egreen_img, mask=mask)\n",
        "plt.imshow(masked_egreen_img)\n",
        "mean_egreen = np.nanmean(egreen_img[mask > 0])\n",
        "print(f\"egreen: {mean_egreen}\")"
      ],
      "metadata": {
        "id": "HGSmAxgMd8OW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Which other indexes you think will be interesting to analyse? Why?\n",
        "To give you a hint, our cropreporter meassures at the following wavelengths:\n",
        "- [RGB wavelengths](https://github.com/danforthcenter/plantcv/blob/10a33cc8cdba9c7a5f8dda0637e4172cbf6681b5/plantcv/plantcv/photosynthesis/read_cropreporter.py#L282C9-L282C21)\n",
        "- [Spectral wavelengths](https://github.com/danforthcenter/plantcv/blob/10a33cc8cdba9c7a5f8dda0637e4172cbf6681b5/plantcv/plantcv/photosynthesis/read_cropreporter.py#L289C14-L289C14)\n"
      ],
      "metadata": {
        "id": "1r4ii0Rq8I1U"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Photosynthesis analysis\n",
        "Unfortunately, there is no consensus on how photosythesis measurement should be done or how to name the frames, as a result different terms are used in literature. For the cropreporter PAM measurements it is as follows:\n",
        "\n",
        "PAM light\n",
        "- **Flight** Correction frame (actinic light on)\n",
        "- **Fp** Fluorescence level light adapted plant of averaged measurement pulses,\n",
        "already corrected for Flightt\n",
        "- **Fmp** Fluorescence level during saturation around 0.8s of averaged measurement pulses,\n",
        "already corrected for Fs\n",
        "- **Fs** Correction frame (actinic light and saturation light on)\n",
        "\n",
        "PAM dark\n",
        "- **Fdark** Correction frame (no light)\n",
        "- **F0** Fluorescence level dark adapted plant of averaged measurement pulses,\n",
        "already corrected for Fdark\n",
        "- **Fm** Fluorescence level during saturation around 0.8s of averaged measurement pulses,\n",
        " already corrected for Fs\n",
        "- **Fs** Correction frame (saturation light on)\n",
        "\n",
        "We can compute the Fp/Fmp or F0/Fm, also known as Fv/Fm and Fq'/Fm'.\n",
        "The raw measurements (Fm+Fs) are not stored and cannot be used for analysis.\n",
        "\n",
        "#### Note, there are strange values in the raw Fm an Fv data of cropreporter, which results in questionable FV/FM values. This is known issue and not solved yet."
      ],
      "metadata": {
        "id": "t3-vhiO2LfHR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# In our meassurements we ofter see zero Fm values for plant pixels.\n",
        "# Therefore we implemented this data quality check: How many pixels in the Fm/Fmp image are zero?\n",
        "FM = None\n",
        "if not ps.pam_dark is None:\n",
        "  #Extract FM data as numpy array\n",
        "  FM = ps.pam_dark.sel(frame_label='Fm').to_numpy()\n",
        "elif not ps.pam_light is None:\n",
        "  #Extract FMp data as numpy array\n",
        "  FM = ps.pam_dark.sel(frame_label='Fmp').to_numpy()\n",
        "elif not ps.ojip_dark is None:\n",
        "  #Extract FM data as numpy array\n",
        "  FM = ps.ojip_dark.sel(frame_label='Fm').to_numpy()\n",
        "\n",
        "if not FM is None:\n",
        "  output = np.full(FM.shape, np.nan)\n",
        "  output[mask > 0] = FM[mask > 0]\n",
        "  output = output.flatten()\n",
        "  print(f\"{len(output[output == 0])} of the {len(output)} pixels, {int(len(output[output == 0])/len(output)*100)}% in the FM image are zero\")"
      ],
      "metadata": {
        "id": "wkh1eD-eLeKk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Calculate and analyze PSII efficiency\n",
        "if not ps.pam_dark is None:\n",
        "  yii_global, yii_chart = pcv.analyze.yii(ps_da=ps.pam_dark, labeled_mask=mask, measurement_labels=[\"Fv/Fm\"])\n",
        "elif not ps.pam_light is None:\n",
        "  yii_global, yii_chart = pcv.analyze.yii(ps_da=ps.pam_light, labeled_mask=mask, measurement_labels=[\"Fq'/Fm'\"])\n",
        "elif not ps.ojip_dark is None:\n",
        "  #Extract FM data as numpy array\n",
        "  yii_global, yii_chart = pcv.analyze.yii(ps_da=ps.ojip_dark, labeled_mask=mask, measurement_labels=[\"Fv/Fm\"])\n",
        "\n",
        "# yii is how plantcv calls the fv/fm or Fq'/Fm'.\n",
        "stats = {'min_yii': round(float(yii_global.min().values),2), 'mean_yii': round(float(yii_global.mean().values),2), 'max_yii': round(float(yii_global.max().values),2)}\n",
        "stats"
      ],
      "metadata": {
        "id": "mojOV0BYmfz9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "qV7iuX1ftVDq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "rVGn-oB4rVPu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Questions\n",
        "- Can you think of more quality checks to perform on the data?\n",
        "- Which other parameters can you compute from this data?"
      ],
      "metadata": {
        "id": "uNCnxKjA-fw_"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "XYyia1uNn9Lv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Example to process an dataset\n",
        "This notebook gives examples aof what can be done with plantcv using a single plant.\n",
        "Once you hav found the parameters you want to analyse you might want to analyse all the data of an entire experiment. In this section an example workflow is given:"
      ],
      "metadata": {
        "id": "THymNuRUL9j8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Preperation\n",
        "from glob import glob\n",
        "from pathlib import Path\n",
        "from json import load\n",
        "import pandas as pd\n",
        "pcv.params.debug = \"None\""
      ],
      "metadata": {
        "id": "QgVg_aXcifC3"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "all_parameters_list = []\n",
        "datafolder2 = \"/content/drive/My Drive/P4P_workshop_data/Lucia\"\n",
        "# Loop over all plants\n",
        "for filename in glob(datafolder2 + \"/**/*.INF\", recursive=True):\n",
        "  # Load data\n",
        "  print(f'analyzing: {filename}')\n",
        "  a_ps = pcv.photosynthesis.read_cropreporter(filename=filename)\n",
        "\n",
        "  # Load metadata\n",
        "  path = Path(filename)\n",
        "  metadata_filename = path.parent.parent.absolute().joinpath('metadata.json')\n",
        "  with open(metadata_filename) as file:\n",
        "      metadata = load(file)\n",
        "\n",
        "  # Mask\n",
        "  mask = pcv.threshold.otsu(gray_img=img_as_ubyte(ps.chlorophyll.sel(frame_label=\"Chl\").data), object_type=\"light\")\n",
        "\n",
        "  # Correct wavelengths\n",
        "  ps.spectral.wavelength_dict = {475: 0.0, 540: 1.0, 550: 2.0, 640: 3.0, 710: 4.0, 770: 5.0}\n",
        "\n",
        "  # NDVI\n",
        "  ndvi_img = pcv.spectral_index.ndvi(hsi=ps.spectral, distance=30).array_data\n",
        "  masked_ndvi_img = cv2.bitwise_and(ndvi_img, ndvi_img, mask=mask)\n",
        "  mean_ndvi = round(np.nanmean(ndvi_img[mask > 0]),2)\n",
        "  # Optional add other traits\n",
        "  #....\n",
        "\n",
        "  #....\n",
        "  # Create dictionary with plant parameters\n",
        "  plant_parameters = {'filename': filename,\n",
        "                      'Datetime': metadata['Datetime'],\n",
        "                      'ExperimentId': metadata['ExperimentId'],\n",
        "                      'TreatmentId': metadata['TreatmentId'],\n",
        "                      'PlantId': metadata['PlantId'],\n",
        "                      'CameraHeight': metadata['Extra']['Height'],\n",
        "                      'ConfigFile': metadata['Extra']['Height'],\n",
        "                      'mean_ndvi': mean_ndvi\n",
        "                      }\n",
        "  all_parameters_list.append(plant_parameters)"
      ],
      "metadata": {
        "id": "-Auo-G9LMBMh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Combine all dictionaries into one table\n",
        "parameter_table = pd.DataFrame(all_parameters_list)\n",
        "parameter_table"
      ],
      "metadata": {
        "id": "-Ows3NyujX5_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Save table as csv\n",
        "parameter_table.to_csv('extracted_cropreporter_traits.csv', index = False)"
      ],
      "metadata": {
        "id": "hguIAY3hj1AF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Bonus: Plotting\n",
        "Comparing data is always nicer with some graphs, while this can be done in Excel. Tools like plotly enable you to make fancy looking interactive graphs.\n",
        "For more information see: https://plotly.com/graphing-libraries/"
      ],
      "metadata": {
        "id": "vl1395aGHzOR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import plotly.graph_objects as go\n",
        "import plotly.express as px\n",
        "import pandas as pd\n",
        "from plotly.subplots import make_subplots"
      ],
      "metadata": {
        "id": "iG1UATxgITAG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Load a prepared dataset\n",
        "croprep_df = pd.read_csv(datafolder + '/cropreporter_traits.csv')\n",
        "# Convert datetime\n",
        "croprep_df['Datetime'] = pd.to_datetime(croprep_df['Datetime'], infer_datetime_format=True)\n",
        "\n",
        "croprep_df.head(5)"
      ],
      "metadata": {
        "id": "_odJf5gRJJZi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "fig = make_subplots()\n",
        "\n",
        "for plantid in croprep_df['PlantId'].unique():\n",
        "  plant_df = croprep_df.loc[croprep_df['PlantId'] == plantid]\n",
        "  fig.add_trace(go.Scatter(x = plant_df['Datetime'], y=plant_df['MeanEgreen'], mode = 'lines', name = plantid))\n",
        "\n",
        "# Set titel\n",
        "fig.update_layout(title_text=\"Mean Egreen of all plants\",\n",
        "                   xaxis_title=\"Datetime\",\n",
        "                  yaxis_title=\"Egreen\",\n",
        "                  legend_title=\"Plants\",)\n",
        "# Interactivity\n",
        "fig.update_layout(\n",
        "    hovermode=\"x unified\", #show trace line\n",
        "    xaxis=dict( # Buttons at the top\n",
        "    rangeselector=dict(\n",
        "        buttons=list([\n",
        "            dict(count=1, label=\"day\", step=\"day\",stepmode=\"backward\"),\n",
        "            dict(step=\"all\")\n",
        "        ])\n",
        "    ),\n",
        "    rangeslider = dict(visible=True), type=\"date\")\n",
        ")\n",
        "fig.show()"
      ],
      "metadata": {
        "id": "xMPMA4t9KeN9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Save the graph to disk\n",
        "fig.write_html('climate plot.html')"
      ],
      "metadata": {
        "id": "QUlDIgtLMX35"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Bonus two: analyse directly on Azure\n",
        "On a regular basis NPEC is contacted by researchers which simply don't have enough diskspace to download their entire dataset.\n",
        "As NPEC we run into the same issue, therefore we run the analysis directly on the data in Azure, without downloading it first!\n",
        "\n",
        "We had to slightly modify plantcv to work in this way, the modified git is here:\n",
        "https://github.com/wurDevTim/plantcv/tree/azure\n",
        "And can be installed by chaning\n",
        "\n",
        "```\n",
        "!pip install plantcv\n",
        "```\n",
        "to:\n",
        "\n",
        "\n",
        "```\n",
        "!pip install git+https://github.com/wurDevTim/plantcv.git@azure\n",
        "```\n",
        "\n"
      ],
      "metadata": {
        "id": "Fauyw2w8PqNc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install azure-storage-blob\n",
        "# In this notebook we already have plantcv, this must be uninstalled first\n",
        "!pip uninstall plantcv\n",
        "!pip install git+https://github.com/wurDevTim/plantcv.git@azure"
      ],
      "metadata": {
        "id": "ujuGydrARHpZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Imports\n",
        "import pandas as pd\n",
        "from json import loads\n",
        "from plantcv import plantcv as pcv\n",
        "from azure.storage.blob import ContainerClient"
      ],
      "metadata": {
        "id": "CNLbW-mtRGdh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "It's not advised to share access keys in a git, not even read only only keys to a test environment like we are using here. Therefore please ask Tim for a test key if want to try this!\n",
        "\n",
        "In this example we use the SAS URL (Shared access token) as used in the Azure storage explorer.\n",
        "\n",
        "The used Permissions are Read & List\n"
      ],
      "metadata": {
        "id": "9DjMLr4EUvPc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "sas_url = \"<paste sas url token here>\"\n",
        "container_client = ContainerClient.from_container_url(sas_url)"
      ],
      "metadata": {
        "id": "LkUIZMqAQ2SL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "This version of pandas needs the exact file locations, therefore we have to do a little magic:"
      ],
      "metadata": {
        "id": "aknGqWkrZ-Rq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Retrieve a list with all cropreporter files in the storage account\n",
        "blob_list = container_client.list_blobs()\n",
        "filenames = []\n",
        "for blob in blob_list:\n",
        "  filename = blob.name\n",
        "  # Only cropreporter data\n",
        "  if filename.split('/',1)[0] == \"CropReporter\":\n",
        "    filenames.append(filename)"
      ],
      "metadata": {
        "id": "TxRnORWHWTWD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Switch to pandas for more advanced sorting options\n",
        "filenames_df = pd.DataFrame({'filename':filenames})\n",
        "#Split into folder & filename:\n",
        "filenames_df[['folder','data','file']] = filenames_df['filename'].str.rsplit(\"/\", n=2, expand=True)\n",
        "\n",
        "# Modify the foldernames of the metadata so the foldernames match the names of the files\n",
        "# The data is in a subfolder.\n",
        "filenames_df.loc[filenames_df['file'] == 'metadata.json', 'folder'] = filenames_df.loc[filenames_df['file'] == 'metadata.json', 'folder'] + '/' + filenames_df.loc[filenames_df['file'] == 'metadata.json', 'data']\n",
        "\n",
        "filenames_df"
      ],
      "metadata": {
        "id": "FwfRSZmrZOKM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Loop over all folders and extract the files:\n",
        "for foldername in filenames_df['folder'].unique():\n",
        "  # One meassurement\n",
        "  kwargs = {}\n",
        "  for idx, row in filenames_df.loc[filenames_df['folder'] == foldername].iterrows():\n",
        "    if row['file'] == 'metadata.json':\n",
        "      metadata_file_path = row['filename']\n",
        "    elif 'HDR' in row['file']:\n",
        "      kwargs['hdr_blob_name'] = row['filename']\n",
        "    elif 'CHL' in row['file']:\n",
        "      kwargs['chl_blob_name'] = row['filename']\n",
        "    elif 'CLR' in row['file']:\n",
        "      kwargs[\"clr_blob_name\"] = row['filename']\n",
        "    elif 'SPC' in row['file']:\n",
        "      kwargs[\"spc_blob_name\"] = row['filename']\n",
        "    elif 'PMD' in row['file']:\n",
        "      kwargs[\"pmd_blob_name\"] = row['filename']\n",
        "    elif 'PML' in row['file']:\n",
        "      kwargs[\"pml_blob_name\"] = row['filename']\n",
        "\n",
        "  # Load the metadata json\n",
        "  blob_client = container_client.get_blob_client(metadata_file_path)\n",
        "  stream = blob_client.download_blob(encoding='UTF-8')\n",
        "  metadata_dict = loads(stream.readall())\n",
        "  # Load plant into plantcv\n",
        "  ps = pcv.photosynthesis.read_cropreporter(container_client, kwargs)\n",
        "  #-----------------------------------\n",
        "  # Process the data, forexample using one of the methods shown above\n",
        "  #-----------------------------------"
      ],
      "metadata": {
        "id": "FyaaPoa1Z7Dk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "F_X33mRGWTuf"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}