<jupyter_start><jupyter_text>1.2. Preprocessing weather indexes<jupyter_code>library("tidyverse")
library("weathercan")
library("lubridate") # yday(), ymd()
library("Amelia")<jupyter_output>Attaching package: 'lubridate'


The following object is masked from 'package:base':

    date


Loading required package: Rcpp

## 
## Amelia II: Multiple Imputation
## (Version 1.7.6, built: 2019-11-24)
## Copyright (C) 2005-2020 <NAME>, <NAME> and <NAME>
## Refer to http://gking.harvard.edu/amelia/ for more information
##<jupyter_text>Load the data frame processed in previous file `1.1_preprocessing-1.ipynb`.<jupyter_code>data_df <- read.csv('output/balanced_potato_df.csv')<jupyter_output><empty_output><jupyter_text>Compute 5 years (mean or total) weather indices Because the process is long, you can load `data_stations.csv` file and make a left join with `data_df` using "Annee", "LatDD", "LonDD", "DatePlantation", and "DateRecolte" as common features.<jupyter_code>data_stations = read_csv("output/weather_stations.csv")<jupyter_output>Parsed with column specification:
cols(
  Annee = [32mcol_double()[39m,
  LatDD = [32mcol_double()[39m,
  LonDD = [32mcol_double()[39m,
  DatePlantation = [34mcol_date(format = "")[39m,
  DateRecolte = [34mcol_date(format = "")[39m,
  station_id = [32mcol_double()[39m,
  temp_moy_5years = [32mcol_double()[39m,
  prec_tot_5years = [32mcol_double()[39m,
  sdi_5years = [32mcol_double()[39m,
  gdd_5years = [32mcol_double()[39m
)<jupyter_text>Custom fonctions<jupyter_code># shannon diversity index
SDI_f <- function(x) { 
  p <- x/sum(x, na.rm = TRUE)
  SDI <- sum(p * log(p), na.rm = TRUE) / log(length(x))
  return(SDI)
}
# growing degree-days
GDD_f <- function(x, delim = 5) {
  sum(x[x >= delim], na.rm = TRUE)
}<jupyter_output><empty_output><jupyter_text>Create a data frame for computations<jupyter_code>data_df$DatePlantation <- ymd(data_df$DatePlantation)
data_df$DateRecolte <- ymd(data_df$DateRecolte)<jupyter_output>Warning message:
" 1 failed to parse."Warning message:
" 1 failed to parse."<jupyter_text>If you've loaded the `weather_stations.csv`, from here, you coud jump to the fith code chaine from the end (insertion into the main table).<jupyter_code>data_stations <- data_df %>%
  distinct(Annee, LatDD, LonDD, DatePlantation, DateRecolte) %>%
  filter(!is.na(Annee))
glimpse(data_stations)<jupyter_output>Observations: 323
Variables: 5
$ Annee          [3m[90m<int>[39m[23m 2003, 2004, 2004, 2004, 2004, 2004, 2004, 2004, 2005...
$ LatDD          [3m[90m<dbl>[39m[23m 46.75306, 46.76444, 48.79722, 48.79722, 48.79722, 48...
$ LonDD          [3m[90m<dbl>[39m[23m -72.33861, -72.33083, -72.24250, -72.24250, -72.2425...
$ DatePlantation [3m[90m<date>[39m[23m 2003-05-19, 2004-05-30, 2004-05-31, 2004-06-05, 200...
$ DateRecolte    [3m[90m<date>[39m[23m 2003-09-17, 2004-09-27, 2004-09-20, 2004-09-15, 200...<jupyter_text>Impute missing dates if year is recorded<jupyter_code># Check missing dates
data_stations %>%
  filter(is.na(DatePlantation))
# Check missing dates
data_stations %>%
  filter(is.na(DateRecolte))
data_stations_imp <- data_stations %>%
                          mutate(DatePlantation_yd = yday(DatePlantation), # yday = year day
                                 DateRecolte_yd = yday(DateRecolte))

data_stations_imp <- amelia(x = data_stations_imp %>%
                         select(Annee, LatDD, LonDD, DatePlantation_yd, 
                                DateRecolte_yd), m=1)$imputations[[1]]
data_stations_imp$Annee <- as.Date(paste0(data_stations_imp$Annee, "-01-01"))

data_stations$DatePlantation <- data_stations_imp$Annee + data_stations_imp$DatePlantation_yd - 1
data_stations$DateRecolte <- data_stations_imp$Annee + data_stations_imp$DateRecolte_yd - 1
data_stations$DatePlantation <- as.Date(data_stations$DatePlantation)
data_stations$DateRecolte <- as.Date(data_stations$DateRecolte)<jupyter_output><empty_output><jupyter_text>Load data from Environnement Canada (weathercan)<jupyter_code>year_step <- 5 # start 5 years preceeding trial year
data_stations$station_id <- NA # identify weather stations
station_weather <- list() # list of weather stations containning climate data

for (i in 1:nrow(data_stations)) {
  print(paste(i, "/", nrow(data_stations), '...'))
  
  # identify the nearest station. Search for stations by name or location
  all_stations <- stations_search(coords = c(data_stations$LatDD[i], 
                                             data_stations$LonDD[i]),  
                                  interval = "day",
                                  dist = 500)
  annee <- data_stations$Annee[i]
  
  # Available station for nearest year
  closest_station <- all_stations %>%
                        filter(start <= annee - year_step) %>%
                        filter(end >= annee) %>%
                        slice(which.min(distance))

  data_stations$station_id[i] <- closest_station$station_id %>% 
                                      as.character() %>% 
                                      as.numeric()
  
  print(paste("Station id:", data_stations$station_id[i]))
  
  # # Download weather data from Environment Canada of the station identified
  station_weather[[i]] <- weather_dl(station_ids = data_stations$station_id[i], 
                                     start = as.Date(paste0(annee - year_step, "-01-01")),
                                     end = as.Date(paste0(annee, "-01-01")),
                                     interval = "day")
  
  print(paste(i, "/", nrow(data_stations)))
}
# name of the list elements
names(station_weather) <- as.character(data_stations$station_id)<jupyter_output>[1] "1 / 323 ..."
[1] "Station id: 5203"
[1] "1 / 323"
[1] "2 / 323 ..."
[1] "Station id: 5203"
[1] "2 / 323"
[1] "3 / 323 ..."
[1] "Station id: 5929"
[1] "3 / 323"
[1] "4 / 323 ..."
[1] "Station id: 5929"
[1] "4 / 323"
[1] "5 / 323 ..."
[1] "Station id: 5929"
[1] "5 / 323"
[1] "6 / 323 ..."
[1] "Station id: 5929"
[1] "6 / 323"
[1] "7 / 323 ..."
[1] "Station id: 5929"
[1] "7 / 323"
[1] "8 / 323 ..."
[1] "Station id: 5929"
[1] "8 / 323"
[1] "9 / 323 ..."
[1] "Station id: 5203"
[1] "9 / 323"
[1] "10 / 323 ..."
[1] "Station id: 5203"
[1] "10 / 323"
[1] "11 / 323 ..."
[1] "Station id: 5222"
[1] "11 / 323"
[1] "12 / 323 ..."
[1] "Station id: 5222"
[1] "12 / 323"
[1] "13 / 323 ..."
[1] "Station id: 5222"
[1] "13 / 323"
[1] "14 / 323 ..."
[1] "Station id: 5929"
[1] "14 / 323"
[1] "15 / 323 ..."
[1] "Station id: 5255"
[1] "15 / 323"
[1] "16 / 323 ..."
[1] "Station id: 5929"
[1] "16 / 323"
[1] "17 / 323 ..."
[1] "Station id: 5255"
[1] "17 / 323"
[1] "18 / 323 ..."
[1] "Station id: 5255"
[1] "18[...]<jupyter_text>Compute indexes for each station<jupyter_code># Initialize columns to record indexes
data_stations$temp_moy_5years <- NA
data_stations$prec_tot_5years <- NA
data_stations$sdi_5years <- NA
data_stations$gdd_5years <- NA

for (i in 1:nrow(data_stations)) {
  starts <- data_stations$DatePlantation[i] - (year_step:1)*365.25 # beginings of seasons
  ends <- data_stations$DateRecolte[i] - (year_step:1)*365.25 # ends of seasons, for previous years
  # Create vectors to record indices from previous years
  temp_moy_j <- rep(NA, length(starts))
  prectot_j <- temp_moy_j
  sdi_j <- temp_moy_j
  gdd_j <- temp_moy_j
  for (j in 1:length(starts)) {
    # Filter the downloaded weather data tables to keep only data measured between the season dates
    season <- station_weather[[i]] %>%
      filter(date >= starts[j]) %>%
      filter(date <= ends[j])
    temp_moy_j[j] <- mean(season$mean_temp, na.rm = TRUE) # average temperature of the season
    prectot_j[j] <- sum(season$total_precip, na.rm = TRUE) # total precipitations of the season
    sdi_j[j] <- SDI_f(season$total_precip) # season SDI
    gdd_j[j] <- GDD_f(season$mean_temp, delim = 5) # season GDD
  }
  data_stations$temp_moy_5years[i] <- mean(temp_moy_j, na.rm = TRUE) # average temperature
  data_stations$prec_tot_5years[i] <- mean(prectot_j, na.rm = TRUE) # average total rainfalls
  data_stations$sdi_5years[i] <- mean(sdi_j, na.rm = TRUE) # average SDI
  data_stations$gdd_5years[i] <- mean(gdd_j, na.rm = TRUE) # average GDD
}
write_csv(x = data_stations, path = "output/weather_stations.csv") # a backup<jupyter_output><empty_output><jupyter_text>If you've loaded the `weather_stations.csv`, continue from here.  Insert into the main table<jupyter_code>weath_col <- c(
            'Annee', 'LatDD', 'LonDD', 'DatePlantation', 'DateRecolte', 
            'temp_moy_5years', 'prec_tot_5years', 'sdi_5years', 'gdd_5years'
            )
data.frame(weath_col, weath_col %in% colnames(data_df))
data_df <- left_join(data_df, data_stations[weath_col], 
                     by = c("Annee", "LatDD", "LonDD", "DatePlantation", "DateRecolte"))
wcol <- c('temp_moy_5years', 'prec_tot_5years', 'sdi_5years', 'gdd_5years')
data.frame(wcol, wcol %in% colnames(data_df))
write_csv(data_df, "output/pr_potato_df.csv") # backup<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15518773605865363, 'fraction_numerical': 0.06987336147522773, 'mean_word_length': 3.5770208439247586, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '503809', 'license': 'MIT', 'n_tokens_mistral': 3973, 'n_tokens_neox': 3370, 'n_words': 748, 'path': '1.2_preprocessing-weather-indexes.ipynb', 'repo_name': 'rgoals/2019_Site-specific-potato-npk-model', 'size': '9023'}