{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.6.0","file_extension":".r","codemirror_mode":"r"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81000,"databundleVersionId":8812083,"sourceType":"competition"}],"dockerImageVersionId":30749,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Context\nI'm only using Kaggle Notebooks to run scripts (and everything I do, I'm making public), so I'll separate the steps in order to save memory and run time. The rationale of the functions and divisions used here is explained in this tutorial: https://www.kaggle.com/code/picsoflily/tutorial-weather-aggregation-using-thermal-sum But the code on this script is an improved version of that one.\n\nThis is an updated version of the R_Script01 that allows for different thresholds simultaneously, i.e., one in which the cycle is divided in three periods (form 0 to 500 to 1000 to 2000) and one in which it's divided in two (from 0 to 1000 to 2000).","metadata":{"_uuid":"fd90dd70-c20a-4f31-ae0e-6a4909286063","_cell_guid":"d855b14d-d921-497d-bb39-e62bfc682973","trusted":true}},{"cell_type":"markdown","source":"# Libraries","metadata":{"_uuid":"be4bb095-0f8b-4780-8c7c-d68730a35239","_cell_guid":"bf1d5d07-4f0b-4d85-b1c0-44d518a0b382","trusted":true}},{"cell_type":"code","source":"options(repr.plot.width = 30, repr.plot.height = 22)\nlibraries <- c(\"here\", \"tidyverse\", \"data.table\", \"zoo\", \"arrow\")\n\nlapply(libraries, require, character.only = TRUE)","metadata":{"_uuid":"b41bfdce-39df-4cbc-8d8d-74231979c0fc","_cell_guid":"4787e687-a79c-4376-938a-c309ee633b9e","collapsed":false,"_kg_hide-output":true,"_kg_hide-input":true,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Functions","metadata":{"_uuid":"036576f5-793b-4441-b5d1-780463c0f59c","_cell_guid":"11395600-4bd0-4876-bbbd-0b51c71b6a83","trusted":true}},{"cell_type":"markdown","source":"## Feature engineering functions","metadata":{"_uuid":"3a18aceb-3454-41bc-bc5c-c3f3c37b2406","_cell_guid":"3c9cafe5-8a38-48e0-8bcb-11cc1bf64c12","trusted":true}},{"cell_type":"code","source":"def_tt <- function(x, tbase, tupper, thresholds){\n    \n    # Returns das for each predefined threshold of accumulated thermal time in the cycle\n    # x: ID and 30 to 209\n    # x_: ID, das (days after simulation started) and correspondent accumulated tt\n    \n    # Does not differentiate between \"never crossed the threshold because temperatures were low\" and \"cycle ended before reaching the following threshold\"\n    \n    # tbase <- 8\n    \n    # thresholds <- c(1000, 2000)\n    thresholds_name <- paste0(\"tt\", thresholds)\n    \n    # Calculates thermal time for each day and accumulated thermal time for each ID\n    x_tt <- data.table::melt(x, id.vars = \"ID\", variable.name = \"das\",\n                             variable.factor = FALSE)\n    x_tt[, das := as.numeric(das)]\n    \n    x_tt[, value := ifelse(value > tupper, 0, value)]\n    x_tt[, value := ifelse(value - tbase < 0, 0, value - tbase)]\n    \n    x_tt <- x_tt[order(ID, das)][, tt := cumsum(value),\n                                 keyby = c(\"ID\")][, value := NULL]\n    x_tt[, (thresholds_name) := as.list(thresholds\n                                        # 0.1*max(tt), 0.25*max(tt), 1/3)*max(tt), 0.5*max(tt), 800\n    ),\n    keyby = ID]\n    \n    # Identifies dates closest to the thresholds\n    x_tt <- data.table::melt(x_tt,\n                             id.vars = c(\"ID\", \"das\", \"tt\"),\n                             measure.vars = thresholds_name,\n                             variable.name = \"accum_tt\")\n    x_tt[, dif:= abs(tt-value)][, value := NULL]\n    x_tt[, flag := (dif == min(dif)), keyby = c(\"ID\", \"accum_tt\")]\n    \n    x_tt <- x_tt[flag == TRUE]\n    x_tt <- x_tt[x_tt[, .I[flag][1], by = .(ID, accum_tt)][, V1]]\n    x_tt <- x_tt[, .(ID, das, tt, accum_tt)]\n    x_tt <- unique(x_tt, by = c(\"das\", \"ID\"))\n    setkey(x_tt, NULL)\n\n    return(x_tt)\n    \n}\n\nfeat_eng <- function(x, tt_das, var_, funcs){\n    \n    # Returns aggregated time-series for each threshold of tt of the variable\n    # x: ID and 30 to 239\n    # tt_das: ID, das and correspondent tt\n    # var_: variable to which dataset refers\n    \n    x_long <- data.table::melt(x, id.vars = \"ID\", variable.name = \"das\",\n                               variable.factor = FALSE)\n    x_long <- x_long[, das := as.numeric(das)][order(ID, das)]\n       \n    if (!is.null(tt_das)){\n        \n        # Create groups by thresholded period\n        \n        x_mod <- merge(x_long, tt_das, all = TRUE)[, tt := NULL]\n        rm(x_long)\n        x_mod[, accum_tt := as.numeric(accum_tt)]\n        # We will ascribe the groups from last day to first, so we use 0 as a\n        # placeholder for group code to the last day if it exceeded the larger threshold,\n        # carry the value until the previous group is found code and\n        # end x_mod with ID, das, tt value and accum_tt group\n        x_mod[, accum_tt := ifelse(das == max(das) & is.na(accum_tt), 0, accum_tt),\n              keyby = \"ID\"]\n        x_mod[, accum_tt := zoo::na.locf(accum_tt, fromLast=TRUE)]\n        x_mod <- x_mod[accum_tt != 0]\n        \n        x_all <- data.table(unique(x_mod[, \"ID\"]))\n        setkey(x_all, NULL)\n        \n        if(var_ == \"temp\"){\n            # We then add, just for analysis, which day of the cycle corresponds\n            # to the date of each threshold\n            tt_das <- data.table::dcast(tt_das, ID ~ accum_tt,\n                                        value.var = \"das\")[, ID := NULL]\n            x_all <- cbind(x_all, tt_das)\n            \n            # We also add indicators to account for missing values\n            flag_period <- unique(x_mod[, .(ID, accum_tt)])[, temp := 1]\n            flag_period <- data.table::dcast(flag_period, ID ~ accum_tt,\n                                             value.var = \"temp\", fill = 0)\n            \n            cols_ <- setdiff(names(flag_period), \"ID\")\n            new_names <- paste0(\"flag_period_\", cols_)\n            setnames(flag_period, old = cols_, new = new_names)\n            \n            x_all <- cbind(x_all, flag_period[, ID := NULL])\n        }\n        \n        rm(tt_das)\n        \n        # Generates aggregation\n        for (func_name in names(funcs)) {\n            # Here we generate the aggregated series for each period and\n            # join them to the aggregation for the whole season\n            x_periods <- x_mod[, .(funcs[[func_name]](value)), by = c(\"ID\", \"accum_tt\")]\n            x_periods[, new_name := paste(var_, func_name,\n                                          accum_tt, sep = \"_\")][, accum_tt := NULL]\n            x_periods <- data.table::dcast(x_periods, \n                                                  ID ~ new_name, value.var = \"V1\")\n            setkey(x_periods, NULL)\n            \n            x_full <- x_mod[, .(funcs[[func_name]](value)), by = c(\"ID\")]\n            setkey(x_full, NULL)\n            x_full[, new_name := paste(var_, func_name, \"season\", sep = \"_\")]\n            x_full <- data.table::dcast(x_full, ID ~ new_name, value.var = \"V1\")\n            setkey(x_full, NULL)\n            \n            temp <- cbind(x_periods, x_full[, ID := NULL])\n            x_all <- merge(x_all, temp, all = TRUE)\n            setkey(x_all, NULL)\n            \n        }\n        \n    } else {\n        \n        # Aggregation not performed using thermal time info (e.g., precipitation before growth)\n        x_mod <- x_long\n        rm(x_long)\n        x_all <- data.table(unique(x_mod[, \"ID\"]))\n        for (func_name in names(funcs)) {\n            x_full <- x_mod[, .(funcs[[func_name]](value)), by = c(\"ID\")]\n            x_full[, new_name := paste0(func_name, \"_season\")]\n            x_full <- data.table::dcast(x_full, ID ~ new_name, value.var = \"V1\")\n            \n            x_all <- merge(x_all, x_full, all = TRUE)\n            setkey(x_all, NULL)\n            \n        }\n        \n    }\n    \n    return(x_all)\n    \n}\n\nagg_funs <- function(type_, variable_){\n    \n    if(type_ == \"simple\"){\n        \n        if(variable_ == \"tas\"){\n            funcs <- list(\"min\" = min, \"max\" = max, \"sd\" = sd)\n        } else if (variable_ == \"tasmax\"){\n            funcs <- list(\"max\" = max)\n        } else if (variable_ == \"tasmin\"){\n            funcs <- list(\"min\" = min)\n        } else if (variable_ == \"rsds\"){\n            funcs <- list(\"sum\" = sum)\n        } else if (variable_ == \"pr\"){\n            funcs <- list(\"sum\" = sum)\n        } else {\n            funcs <- NULL\n        }\n        \n    } else {\n        \n        if(variable_ == \"tas\"){\n            funcs <- list(\"avg\" = mean, \"min\" = min, \"max\" = max, \"sd\" = sd)\n        } else if (variable_ == \"tasmax\"){\n            funcs <- list(\"avg\" = mean, \"max\" = max, \"sd\" = sd)\n        } else if (variable_ == \"tasmin\"){\n            funcs <- list(\"avg\" = mean, \"min\" = min, \"sd\" = sd)\n        } else if (variable_ == \"rsds\"){\n            funcs <- list(\"avg\" = mean, \"sum\" = sum)\n        } else if (variable_ == \"pr\"){\n            funcs <- list(\"sum\" = sum, \"sd\" = sd)\n        } else {\n            funcs <- NULL\n        }\n        \n    }\n    \n    return(funcs)\n}","metadata":{"_uuid":"178e4939-207b-4461-b793-a5dce6b23d78","_cell_guid":"e24860d5-4c40-4ccd-bc70-121ae2d01c79","collapsed":false,"jupyter":{"outputs_hidden":false},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reading and organizing data functions","metadata":{"_uuid":"c055d4a0-dd1d-43a4-9967-97e278c8876c","_cell_guid":"efe4e496-9cd6-4491-86e8-2beafd56249f","trusted":true}},{"cell_type":"code","source":"read_data <- function(filename, files_, max_rows, cols_sel=NULL, ids=NULL){\n    \n    if(is.null(cols_sel)) {\n        \n        if(is.null(max_rows)){\n            \n            x <- read_parquet(files_[grepl(filename, files_)]) %>%\n                as.data.table()\n        } else {\n            x <- open_dataset(files_[grepl(filename, files_)]) %>%\n                head(n=max_rows) %>%\n                compute() %>%\n                as.data.table()\n        }\n        \n    } else {\n        \n        if(is.null(max_rows)){\n            x <- read_parquet(files_[grepl(filename, files_)],\n                              col_select=cols_sel) %>%\n                as.data.table()\n        } else {\n            x <- open_dataset(files_[grepl(filename, files_)]) %>%\n                head(n=max_rows) %>%\n                compute() %>%\n                as.data.table()\n            \n            x <- x[, ..cols_sel]\n        }\n        \n    }\n    \n    if(!is.null(ids)){\n        x <- x[ID %in% ids]\n    }\n    \n    return(x)\n    \n}\n\ngenerate_set <- function(crop, set, variable, funcs, max_rows, files_,\n                         tt_das=NULL, tbase=8, tupper=40, thresholds=c(\"1000\", \"2000\")){\n    \n    list_out <- list()\n    \n    cols_growth <- c(\"ID\", paste0(c(as.character(seq(30,239)))))\n    cols_pre <- c(\"ID\", paste0(c(as.character(seq(1,30)-1))))\n    \n    varname <- case_when(variable == \"tas\" ~ \"temp\",\n                         variable == \"tasmax\" ~ \"tmax\",\n                         variable == \"tasmin\" ~ \"tmin\",\n                         variable == \"rsds\" ~ \"srad\",\n                         variable == \"pr\" ~ \"prec\")\n    \n    filename <- paste(variable, crop, set, sep = \"_\")\n    \n    if(grepl(\"soil\", variable)){\n        cols_sel_ <- NULL\n    } else if (grepl(\"pr\", variable) & (is.null(tt_das))){\n        cols_sel_ <- cols_pre \n    } else {\n        cols_sel_ <- cols_growth\n    }\n    \n    if (variable != \"solutions\"){\n        \n        x <- read_data(filename, files_, max_rows, cols_sel=cols_sel_)\n        \n        # Read files and processes, in batches, determining the thermal time thresholds\n        if (grepl(\"tas\", variable) & (is.null(tt_das))){\n            # Thermal time\n            tt_das_l <- list()\n            steps_it <- split(x$ID, cut(seq_along(x$ID), 10, labels = FALSE))\n            \n            rm(x)\n            \n            for (it in 1:length(steps_it)){\n                \n                ids <- steps_it[[it]]\n                x_it <- read_data(filename, files_, max_rows, cols_sel=cols_sel_, ids=ids)\n                \n                tt_das_l[[it]] <- def_tt(x_it, tbase, tupper, thresholds)\n                rm(x_it)\n                gc()\n                \n            }\n            \n            list_out[[\"tt_das\"]] <- rbindlist(tt_das_l)\n            \n            return(list_out)\n        }\n        \n        if(grepl(\"soil\", variable)){\n            \n            x_mod <- rename(x, soil=texture_class)\n            \n            # Also batch-read files to generate features\n        } else {\n            \n            x_mod_l <- list()           \n            steps_it <- split(x$ID, cut(seq_along(x$ID), 10, labels = FALSE))\n            \n            rm(x)\n            \n            for (it in 1:length(steps_it)){\n                \n                ids <- steps_it[[it]]\n                x_it <- read_data(filename, files_, max_rows, cols_sel=cols_sel_, ids)\n                \n                tt_das_it <- tt_das\n                if(!is.null(tt_das)){\n                    tt_das_it <- tt_das[ID %in% ids]\n                }\n                \n                x_mod_l[[it]] <- feat_eng(x_it, tt_das_it, varname, funcs)\n                rm(x_it)\n                gc()\n                \n                #cat(paste0(\"it \", it, \" in \", length(steps_it)))\n                \n            }\n            \n            x_mod <- rbindlist(x_mod_l, fill=TRUE)\n            rm(x_mod_l)\n            \n            if(any(grepl(\"tt[0-9]+\", colnames(x_mod)))){\n                var_rem <- grep(\"tt[0-9]+\", colnames(x_mod), value = TRUE)\n                x_mod <- x_mod[, !..var_rem]\n            }\n        }\n        \n        if(grepl(\"pr\", variable) & (is.null(tt_das))){\n            setnames(x_mod, old = \"sum_season\", new = \"sum_pre\")\n            setnames(x_mod, old = \"sd_season\", new = \"sd_pre\")\n        }\n        \n        list_out[[varname]] <- x_mod\n        \n    } else {\n        \n        if(set == \"train\"){\n            \n        filename <- paste(set, \"solutions\", crop, sep = \"_\")\n        list_out[[varname]] <- read_data(filename, files_, max_rows, cols_sel=NULL)\n            \n        }        \n        \n    }\n    \n    return(list_out)\n    \n}","metadata":{"_uuid":"370512ee-ad22-429d-b104-db8f7b649d43","_cell_guid":"0c444a8e-c34d-4604-a0cd-520511e59a0f","collapsed":false,"jupyter":{"outputs_hidden":false},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Constants","metadata":{"_uuid":"120bcb96-3e4f-4c83-be34-15c80d31b3f4","_cell_guid":"c0e09679-070d-4c60-bd8e-7f752441fbf1","trusted":true}},{"cell_type":"code","source":"max_rows_ <- NULL # Allows for testing a subset of the dataset; NULL reads the whole file\npath_inputs <- \"../input/the-future-crop-challenge/\"\nfiles_ <- list.files(path_inputs, full.names = TRUE, pattern = \".parquet\")","metadata":{"_uuid":"c0ff2981-562e-4336-8da0-1186bc56b56f","_cell_guid":"8c4a84e0-cf60-4a28-82b0-ecbfbe1227a5","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load and process","metadata":{"_uuid":"dac09d81-b0aa-409d-8ab7-a27073c683ec","_cell_guid":"94c84bce-d56a-432d-9ca6-93b16c6cd18f","trusted":true}},{"cell_type":"code","source":"# Two divisions\ntt_thresholds_maize_l <- list(c(300, 750, 1000, 1500, 2000), c(1000, 2050))\ntt_thresholds_wheat_l <- list(c(200, 500, 1000, 1500, 2200), c(500, 2050))\n\nfiles_config <- expand.grid(set = c(\"train\", \"test\"), crop = c(\"maize\", \"wheat\"), stringsAsFactors=FALSE) %>%\n    mutate(filename = if_else(set == \"train\", \"tr\", \"tst\"),\n           filename = paste(filename, crop, sep = \"_\"),\n           n_divs = if_else(crop == \"maize\", length(tt_thresholds_maize_l), length(tt_thresholds_wheat_l)),\n           tbase = if_else(crop == \"maize\", 8, 0),\n           tupper = if_else(crop == \"maize\", 50, 45))\n\n# Must have the same lenght as the threshold list lenght\nfunctions_l <- list(\"full\", \"simple\")\n\nlist_vars_ <- c(\"soil_co2\", \"pr\", \"tas\", \"tas\", \"tasmax\", \"tasmin\", \"rsds\", \"pr\", \"solutions\")\n# Must have the same lenght as the threshold list lenght\nlist_vars_l <- list(list_vars_, c(\"tas\", \"tas\", \"tasmax\", \"tasmin\", \"pr\"))\n\n# Must have the same lenght as the threshold list lenght\ninclude_season_l <- list(TRUE, FALSE)\n\n# Must have the same lenght as the threshold list lenght\nsuffix <- c(\"\", \"_v2\")","metadata":{"_uuid":"1e48ff8d-f702-4d9f-9b43-26bdf7cf62c0","_cell_guid":"1fd8eee0-5307-4820-824a-f62fc39acf0b","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # One division\n# tt_thresholds_maize_l <- list(c(300, 750, 1000, 1500, 2050))\n# tt_thresholds_wheat_l <- list(c(200, 500, 1000, 1500, 2200))\n\n# files_config <- expand.grid(set = c(\"train\", \"test\"), crop = c(\"maize\", \"wheat\"), stringsAsFactors=FALSE) %>%\n#     mutate(filename = if_else(set == \"train\", \"tr\", \"tst\"),\n#            filename = paste(filename, crop, sep = \"_\"),\n#            n_divs = if_else(crop == \"maize\", length(tt_thresholds_maize_l), length(tt_thresholds_wheat_l)),\n#            tbase = if_else(crop == \"maize\", 8, 0),\n#            tupper = if_else(crop == \"maize\", 50, 45))\n\n# functions_l <- list(\"full\")\n\n# list_vars_ <- c(\"soil_co2\", \"pr\", \"tas\", \"tas\", \"tasmax\", \"tasmin\", \"rsds\", \"pr\", \"solutions\")\n# list_vars_l <- list(list_vars_)\n\n# include_season_l <- list(TRUE)\n\n# suffix <- c(\"\")","metadata":{"_uuid":"78c1b1e8-2cb6-4c57-975e-6595fc7399c3","_cell_guid":"912c542a-2571-4e1a-9e10-5d27102f36a1","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# files_config <- files_config[1, ]\n\n# files_it <- 1\n# thresh_it <- 1","metadata":{"_uuid":"a47ae714-0464-4cbd-a10e-8a034038e288","_cell_guid":"6f6d398b-5590-4951-be51-6658923deec4","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (files_it in 1:nrow(files_config)){\n    \n#     print(files_it)\n    \n    filename_ <- files_config$filename[files_it]\n    crop_ <- files_config$crop[files_it]\n    set_ <- files_config$set[files_it]\n    tbase_ <- files_config$tbase[files_it]\n    tupper_ <- files_config$tupper[files_it]\n    \n    tt_thresholds_l <- tt_thresholds_wheat_l\n    if(crop_ == \"maize\"){\n        tt_thresholds_l <- tt_thresholds_maize_l\n    }\n    \n    tr_file <- NULL\n    \n    for(thresh_it in 1:files_config$n_divs[files_it]){        \n        \n        tt_das_ <- NULL\n        \n        tt_thresholds <- tt_thresholds_l[[thresh_it]]\n        include_season <- include_season_l[[thresh_it]]\n        list_vars <- list_vars_l[[thresh_it]]\n        funcs_type <- functions_l[[thresh_it]]\n        \n        for(variable_ in list_vars){\n            \n            funcs <- agg_funs(type_=funcs_type, variable_)\n            \n            temp <- generate_set(crop=crop_, set=set_, variable=variable_, funcs,\n                                 max_rows=max_rows_, files_, tt_das=tt_das_, tbase=tbase_, tupper=tupper_,\n                                 thresholds=tt_thresholds)\n            \n            if(!is.null(temp[[\"tt_das\"]])){\n                \n                tt_das_ <- temp[[\"tt_das\"]]\n                \n            } else {\n                \n                if(length(temp) > 0){\n                    \n                    temp <- temp[[1]]\n                    \n                    if(!grepl(\"soil\", variable_)){\n                        temp <- temp[, !c(\"ID\")]\n                        setnames(temp, old = colnames(temp), \n                                 new = paste0(colnames(temp), suffix[thresh_it]))\n                    }\n\n                    if(!include_season){\n                        cols_filt <- colnames(temp)[!grepl(\"season\", colnames(temp))]\n                        temp <- temp[, ..cols_filt]\n                    }\n\n                    tr_file <- cbind(tr_file, temp)\n                \n                }         \n            }\n        }    \n        \n    }\n    \n#     print(sort(colnames(tr_file)))\n    gc()\n    save(tr_file, file=here(paste0(filename_, \".RData\")))\n    save(tt_das_, file=here(paste0(\"tt_das_\", filename_,\".RData\")))\n    \n    write.csv(tr_file, file=here(paste0(filename_, \".csv\")), row.names = FALSE, quote = FALSE)\n    \n}","metadata":{"_uuid":"09a97b54-a55b-41b1-9396-5b49f682d6ce","_cell_guid":"ba0bf982-7576-40a9-b286-dfe44986c8d3","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizations","metadata":{"_uuid":"57ffec12-12b2-4d39-bf4c-e8152f2da272","_cell_guid":"fe7bc699-0834-48e5-b5bf-b719b8afaf47","trusted":true}},{"cell_type":"code","source":"files_l <- list.files(\"./\", pattern=\"(tt_das)_.*_(maize|wheat)(.RData)\", full.names=TRUE)\n\nfor(file_it in files_l){\n    load(file_it)\n    \n    max_tt <- tt_das_ %>%\n        group_by(ID) %>%\n        filter(tt == max(tt)) %>%\n        ungroup() %>%\n        select(ID, tt) %>%\n        rename(max_tt = tt)\n    \n    max_das <- tt_das_ %>%\n        group_by(ID) %>%\n        filter(das == max(das)) %>%\n        ungroup() %>%\n        select(ID, das) %>%\n        rename(max_das = das)\n\n    tt_das_ <- tt_das_ %>%\n        select(-tt) %>%\n        pivot_wider(names_from = accum_tt, values_from = das) %>%\n        left_join(max_tt) %>%\n        left_join(max_das)\n    \n    var_name <- paste0(\"tt\",gsub(\".*tt\", \"\", file_it)) %>%    \n        gsub(\".RData\", \"\", .)\n    assign(var_name, tt_das_)    \n}\n\nrm(tt_das_, var_name, max_tt, max_das)","metadata":{"_uuid":"0469317f-4765-41bb-b7a7-725ea9855b80","_cell_guid":"5e3391e0-f744-4fd4-927b-940bd9687c4e","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-08-10T21:55:16.024960Z","iopub.execute_input":"2024-08-10T21:55:16.026650Z","iopub.status.idle":"2024-08-10T21:55:56.488377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Attention: loading .RData files changed!","metadata":{"_uuid":"06f6a7e8-ef8e-4289-8580-ef1c87ee1001","_cell_guid":"46dd09c2-e94b-431a-a3dc-89ab01843b76","trusted":true}},{"cell_type":"code","source":"files_l <- list.files(\"./\", pattern=\"(tr|tst)_(maize|wheat)(.RData)\", full.names=TRUE)\nfor(file_it in files_l){\n   load(file_it)\n    if(!grepl(\"tt_das\", file_it)){\n        var_name <- gsub(\".RData\", \"\", gsub(\".//\", \"\", file_it))\n        assign(var_name, tr_file)\n    }    \n}","metadata":{"_uuid":"6d47496e-724a-49fc-96bb-c6fa36ea905c","_cell_guid":"88d577fa-03a0-4bce-9d2d-de54a7da3573","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-08-10T21:56:23.554313Z","iopub.execute_input":"2024-08-10T21:56:23.556010Z","iopub.status.idle":"2024-08-10T21:56:45.771971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cols_sel <- c(\"ID\", \"lat\", \"lon\", \"real_year\", \"yield\")\n\nmaize_tr <- tt_das_tr_maize %>%\n    left_join(select(tr_maize, any_of(cols_sel)))\n\nwheat_tr <- tt_das_tr_wheat %>%\n    left_join(select(tr_wheat, any_of(cols_sel)))\n\nmaize_tst <- tt_das_tst_maize %>%\n    left_join(select(tst_maize, any_of(cols_sel)))\n\nwheat_tst <- tt_das_tst_wheat %>%\n    left_join(select(tst_wheat, any_of(cols_sel)))","metadata":{"_uuid":"d5d9c7f9-362d-470b-b113-0cfcb8b76e53","_cell_guid":"0a56fa28-f5c3-4a16-af7e-531e8114e3a2","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-08-10T21:56:45.776409Z","iopub.execute_input":"2024-08-10T21:56:45.777930Z","iopub.status.idle":"2024-08-10T21:56:46.392318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Yield","metadata":{"_uuid":"a0a9ea0f-e14d-4d5e-ba18-2201c7321264","_cell_guid":"0bf81df5-28f9-46f3-9a8c-60feb0f3ba90","trusted":true}},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=maize_tr, aes(x=lon, y=lat, color=yield)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"dcecb2d2-0221-4230-ad9f-1eb530bf2039","_cell_guid":"5b07ee1d-f6a3-4836-9b2b-62cfff13fead","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-08-10T21:57:37.239736Z","iopub.execute_input":"2024-08-10T21:57:37.242279Z","iopub.status.idle":"2024-08-10T21:57:59.451514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=wheat_tr, aes(x=lon, y=lat, color=yield)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"0161f33e-1042-468e-9d40-c2f079d935d0","_cell_guid":"e514a1a4-7ce3-4637-9b66-a7251bc1a34d","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-08-10T21:57:59.456105Z","iopub.execute_input":"2024-08-10T21:57:59.458561Z","iopub.status.idle":"2024-08-10T21:58:17.706689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Maximum accumulated thermal time","metadata":{"_uuid":"b7c9a018-674f-43c9-88b4-909f378367c4","_cell_guid":"a020af59-8b56-4dbc-aea0-101fc1142c8b","trusted":true}},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=maize_tr, aes(x=lon, y=lat, color=max_tt)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"83f8c08a-877b-47c9-94ea-f3c042049229","_cell_guid":"b4fa50eb-7c9f-465d-b3fe-88b4837ce9be","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-08-10T21:53:09.323025Z","iopub.execute_input":"2024-08-10T21:53:09.324492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=wheat_tr, aes(x=lon, y=lat, color=max_tt)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"81548f7d-1026-4969-9932-f197c88de165","_cell_guid":"a18a7664-a505-4085-b498-15e9530e09a9","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=maize_tst, aes(x=lon, y=lat, color=max_tt)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"17c5a842-de8e-4df5-9f33-38e704ade2d2","_cell_guid":"d35cce8c-988d-4ca2-aba5-3f2df2ee1243","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=wheat_tst, aes(x=lon, y=lat, color=max_tt)) +\n    scale_color_continuous(type=\"viridis\") +\n    theme_bw()","metadata":{"_uuid":"2f1c2ab0-033c-4b9c-bae0-71fa7ba4cf5b","_cell_guid":"ff032564-3681-4ded-b40a-59d8cd77c3c7","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Season length","metadata":{"_uuid":"b11df228-e45c-4261-99dd-824707a908fd","_cell_guid":"50bfd883-3ef0-4f1c-acf3-de412d7e4e7e","trusted":true}},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=maize_tr, aes(x=lon, y=lat, color=max_das)) +\n    scale_color_continuous(type=\"viridis\", limits = c(50, 240)) +\n    theme_bw()","metadata":{"_uuid":"dbb3c980-1f5e-4254-813a-1858882fc5d5","_cell_guid":"0a994260-299e-4c7d-adeb-06f5785c31df","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=wheat_tr, aes(x=lon, y=lat, color=max_das)) +\n    scale_color_continuous(type=\"viridis\", limits = c(0, 240)) +\n    theme_bw()","metadata":{"_uuid":"5fa27113-bd27-47e5-9c2e-c6cb3079a6e1","_cell_guid":"b98d48c5-c93f-4c0a-8ef5-b82b7daf5df4","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=maize_tst, aes(x=lon, y=lat, color=max_das)) +\n    scale_color_continuous(type=\"viridis\", limits = c(50, 240)) +\n    theme_bw()","metadata":{"_uuid":"7c7d4a40-c3bf-45b6-b843-e44289645cbf","_cell_guid":"4aa72055-3cc7-43e2-a239-b147cf786df1","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot() +\n    facet_wrap(\"real_year\") +\n    geom_point(data=wheat_tst, aes(x=lon, y=lat, color=max_das)) +\n    scale_color_continuous(type=\"viridis\", limits = c(0, 240)) +\n    theme_bw()","metadata":{"_uuid":"eed2d2c8-d8a2-4bcf-a1cf-088eaefc80f4","_cell_guid":"d502e9ca-105b-445c-8a67-6403ce73f8a8","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]}]}