correcting_meassurement_time.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyPRklf6gguAyH1r2Ogu9Iww",
      "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/correcting_meassurement_time.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Correcting for meassurement time\n",
        "Systems like the cropreporter can only meassure one plant at a time, including dark/light adaptation there can be several hours between the meassurement of the first and last plant.\n",
        "\n",
        "As you can imaging a lot can happen in 3 hours: plants can grow, diseases can spread. To correct for this it's advised to use a simple model to 'shift' all meassuring moment to the same time.\n",
        "\n",
        "The biometric group at Wageningen developend an R package to do this and apply other corrections, for example for the location in the greenhouse: https://biometris.github.io/LMMsolver/\n",
        "\n",
        "for more information on splines please see: https://en.wikipedia.org/wiki/Spline_(mathematics)\n",
        "\n",
        "An important note is that the results can differ between R version, we also encountered this when preparing the workshop.\n",
        "\n"
      ],
      "metadata": {
        "id": "u97oKWmG-i53"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## setup\n",
        "Next to python, R code can also be used in colab notebooks.\n",
        "Here we call R from python instead with rpy2 to enable us to use python code as well.\n",
        "An R cell is marked with '%%R' at the top of the cell."
      ],
      "metadata": {
        "id": "WgiH2fK2-rAv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%load_ext rpy2.ipython"
      ],
      "metadata": {
        "id": "L0EtSYXuiLxu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# indicate that you're running R code\n",
        "%%R\n",
        "\n",
        "# Install packages\n",
        "install.packages(\"LMMsolver\")\n"
      ],
      "metadata": {
        "id": "sYtNu_JT-rI9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Import\n",
        "%%R\n",
        "library(LMMsolver)"
      ],
      "metadata": {
        "id": "sqfYkxKY_oqB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Mount google drive - not found an R alternative, using python instead.\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": "6phkogzOcIhN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [],
      "metadata": {
        "id": "ad-yPPVrB-j1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Function which uses the LMM solver to compute the spline.\n",
        "%%R\n",
        "##------------------------------------------------------------\n",
        "# Inputs:\n",
        "# - Dataframe with all the data\n",
        "# - The unique identified for a specific plant\n",
        "# - A list of traits to fir a spline for\n",
        "# Return:\n",
        "# - A list with predictions for all days in the meassuring period\n",
        "##------------------------------------------------------------\n",
        "compute_spline <- function(df, plant_identifier, trait_list)\n",
        "{\n",
        "  ### Fit 1D spline per plant\n",
        "  for (i in c(1:length(unique(df[[plant_identifier]])))){\n",
        "    plant_id = unique(df[[plant_identifier]])[i]\n",
        "    one_plant <- df[df[[plant_identifier]] == plant_id,]\n",
        "    datenum = one_plant[['datenum']]\n",
        "    preddates <- data.frame(datenum = min(one_plant$date):max(one_plant$date))\n",
        "    # Each day has 24*60*60 = 86400 hours\n",
        "    preddates <- preddates * 86400\n",
        "\n",
        "    # Fit 1D spline per trait\n",
        "    for (trait in trait_list){\n",
        "      # Check for inf values\n",
        "      if (sum(is.infinite(one_plant[[trait]])) > 0) {\n",
        "        print(paste('Warning: infinite value encoutered for plant: ', plant_id, ', trait: ',trait))\n",
        "      }\n",
        "      trait_df <- one_plant[!is.infinite(one_plant[[trait]]),]\n",
        "      # Need at least 2 not inf values\n",
        "      # Need at least 2 not inf values\n",
        "      if (nrow(trait_df) > 2) {\n",
        "        # Nan values will be removed, but they do cause warnings.\n",
        "        m1 <- LMMsolve(fixed = as.formula(paste(trait, \"~\", 1)),\n",
        "                       spline = ~spl1D(x = datenum, nseg = 20),\n",
        "                       data = trait_df)\n",
        "        #summary(m1)\n",
        "\n",
        "        # Note, in some cases you might only want to predict from the first meassurement to the last.\n",
        "        # Especially when the first/last one was a nan.\n",
        "        prediction <- obtainSmoothTrend(m1, newdata = preddates,\n",
        "                                        includeIntercept = T)\n",
        "        # Rename ypred column\n",
        "        names(prediction)[names(prediction) == 'ypred'] <- trait\n",
        "      } else {\n",
        "        print(paste('Warning: not enough values to process plant: ', plant_id, ', trait: ',trait))\n",
        "        prediction <- preddates\n",
        "        prediction[trait] <- NA\n",
        "      }\n",
        "      prediction <- prediction[,c(\"datenum\",trait)]\n",
        "      # Combine results\n",
        "      if (trait == trait_list[1]){\n",
        "        plant_predictions <- prediction\n",
        "      } else {\n",
        "        plant_predictions <- merge(plant_predictions, prediction, by='datenum')\n",
        "      }\n",
        "    }\n",
        "    plant_predictions[[plant_identifier]] = plant_id\n",
        "    if (i == 1){\n",
        "      all_predictions <- plant_predictions\n",
        "    } else {\n",
        "      all_predictions <- rbind(all_predictions, plant_predictions)\n",
        "    }\n",
        "  }\n",
        "  return(all_predictions)\n",
        "}"
      ],
      "metadata": {
        "id": "-EAo2gnyAPDz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Example\n",
        "In this example the cropreporter data from Lucia is used, which has been analysed beforehand.\n",
        "\n",
        "Note: Systems like the cropreporter in NPEC use the local time. If your unlucky your experiment contains both winter & summer time. In this case we would advice to switch to UTC."
      ],
      "metadata": {
        "id": "QDzgQ2BkJfBm"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Load the data\n",
        "%%R\n",
        "df <- read.csv('/content/drive/My Drive/P4P_workshop_data/cropreporter_traits.csv', sep = \",\")\n",
        "head(df)"
      ],
      "metadata": {
        "id": "EXucnEuVJc5S"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# Converting to datetime object\n",
        "df[['Datetime']] <- as.POSIXct(df[['Datetime']], format = \"%Y-%m-%d, %H:%M:%OS\", tz=\"Europe/Paris\")\n",
        "\n",
        "\n",
        "# Datenum stored as integer, exact datetime of measurement\n",
        "df[['datenum']] <- as.integer(df[['Datetime']])\n",
        "\n",
        "# Date is multiplied with 86400 to get value at 00:00:00 of each day\n",
        "df[['date']] <- as.numeric(as.Date(df[['Datetime']]))"
      ],
      "metadata": {
        "id": "Ab85ZOMijL5z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "head(df)"
      ],
      "metadata": {
        "id": "ILYddD_9lT1y"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# All columns\n",
        "colnames(df)"
      ],
      "metadata": {
        "id": "shPdnIUckamb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# List of the columns to interpolate\n",
        "trait_list <- list('mean_yii', 'MeanChlorophyll', 'MeanNdvi', 'MeanEgreen', 'MeanPsri', 'MeanAri', 'MeanMari')\n"
      ],
      "metadata": {
        "id": "D_iUrzvykKWb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# Compute the spline\n",
        "predictions = compute_spline(df, 'PlantId', trait_list)"
      ],
      "metadata": {
        "id": "M_hWaoDMkoxn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# Convert time back to datetime, trick question is this datetime now UTC or Amsterdan local time?\n",
        "predictions[['Datetime']] <- as.POSIXct(predictions[['datenum']]\t, origin=\"1970-01-01\", tz = \"UTC\")"
      ],
      "metadata": {
        "id": "UeoIKyDlkr9J"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# Look at some results\n",
        "head(predictions[predictions[['PlantId']] == unique(predictions[['PlantId']])[1],])"
      ],
      "metadata": {
        "id": "mCSV0mNg3CUe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "# Save predictions\n",
        "write.csv(predictions, \"/content/drive/My Drive/P4P_workshop_data/pspline_predictions.csv\", row.names=FALSE)"
      ],
      "metadata": {
        "id": "BHvMphpdktOB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "r64i9YaXl5Qd"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}