From 7731ffa5a37851d4fbcc909eea3b8cddd90f853e Mon Sep 17 00:00:00 2001 From: Henrique Goulart <45360568+dumontgoulart@users.noreply.github.com> Date: Fri, 17 Mar 2023 20:28:42 +0100 Subject: [PATCH] Updating general files after review --- code/create_raster_harvest_area_br.r | 143 ++ code/create_raster_soy_yield_br.R | 195 ++ code/create_raster_soy_yield_us.R | 129 + code/ga_sim_plots.py | 28 +- code/globiom_grid_fun_am.R | 323 +++ code/globiom_harvest_area_converter.R | 42 + code/globiom_mask_generator.R | 50 + code/gtsm_tests.py | 107 + code/ml_predictions_auto_am.py | 419 --- code/paper2_results_comparison_production.py | 4 +- code/paper_2_americas_sim2.py | 109 +- code/soy_harvest_spam_globiom.csv | 2381 ++++++++++++++++++ code/tiff_to_nc.py | 592 +++++ code/wp3_sim.py | 2 +- 14 files changed, 4062 insertions(+), 462 deletions(-) create mode 100644 code/create_raster_harvest_area_br.r create mode 100644 code/create_raster_soy_yield_br.R create mode 100644 code/create_raster_soy_yield_us.R create mode 100644 code/globiom_grid_fun_am.R create mode 100644 code/globiom_harvest_area_converter.R create mode 100644 code/globiom_mask_generator.R create mode 100644 code/gtsm_tests.py create mode 100644 code/soy_harvest_spam_globiom.csv create mode 100644 code/tiff_to_nc.py diff --git a/code/create_raster_harvest_area_br.r b/code/create_raster_harvest_area_br.r new file mode 100644 index 0000000..cf2ff4d --- /dev/null +++ b/code/create_raster_harvest_area_br.r @@ -0,0 +1,143 @@ +library(sidrar) +library(terra) +library(tidyr) +library(dplyr) +library(sf) +library(ncdf4) +library(raster) +library(rgdal) +library(stringr) +setwd("~/PhD/paper_hybrid_agri/data") + +# period.list <- list(1975:1978) +# +# # 112 for yield, 109 for planting area, 216 for harvesting area +# mun.soy.has <- lapply(period.list, function(x) +# get_sidra(1612, variable = 216, +# period = as.character(x), +# geo = "City", +# classific = c("c81"), +# category = list(2713))) +# +# test_2 <- do.call(rbind,mun.soy.has) +# write.csv(test_2,"soy_harv_area_munic_br_1975_1978.csv") + +# mun.soy.has_2 <- do.call(rbind,mun.soy.has) %>% +# rename(ADM2_PCODE = `Município (Código)`) %>% +# dplyr::select(Ano,ADM2_PCODE,Valor) + +mun.soy.has_2 <- read.csv("soy_harv_area_munic_br_1975:1978.csv") +mun.soy.has_2 <- mun.soy.has_2 %>% rename(ADM2_PCODE = 'Município..Código.') %>% + dplyr::select(Ano,ADM2_PCODE,Valor) + +mun.soy.has_2 <- read.csv("soy_harv_area_munic_br_1975_1978.csv") + +mun.soy.has_3 <- read.csv("soy_harv_area_munic_br.csv") +mun.soy.has <- rbind(mun.soy.has_2, mun.soy.has_3) + +mun.soy.has <- mun.soy.has %>% rename(ADM2_PCODE = 'Município..Código.') %>% + dplyr::select(Ano,ADM2_PCODE,Valor) +mun.soy.has <- dplyr::arrange(mun.soy.has , ADM2_PCODE, Ano) + +BRmun <- st_read("GIS/bra_admbnda_adm2_ibge_2020.shp") %>% + mutate(ADM2_PCODE = substr(ADM2_PCODE,3,9)) %>% + mutate_at(c('ADM2_PCODE'), as.numeric) +st_crs(BRmun) <- 4326 +BRmun <- st_transform(BRmun, crs = 4326) + +shp.soy.has <- left_join(BRmun, mun.soy.has) %>% drop_na() + + +# --------------------------------------------------- +# Save as shapefile +# --------------------------------------------------- +# st_write(shp.soy.has, paste0(getwd(), "/", "soy_ibge_br_harvest_area.shp")) + +# Calculate harvest area % +areas <- read.csv("GIS/areas.csv",sep = ";") %>% + rename(ADM2_PCODE = CD_GCMUN) %>% rename(area = AR_MUN_2019) %>% + mutate(area = 100*as.integer(str_replace(area,",","."))) + +area_municipality <- inner_join(areas,shp.soy.has) +area_municipality$harv_area_frac <- area_municipality$Valor / area_municipality$area +area_municipality$Valor[area_municipality$Valor==0] <- NA +area_municipality = area_municipality %>% drop_na() + +# Check for values above 1 - error +#sum(area_municipality$harv_area_frac > 1, na.rm=T) # SHOULD BE ZERO +area_municipality %>% filter(harv_area_frac < 0.01) +sum(area_municipality$harv_area_frac < 0.01, na.rm=T) +#unique(area_municipality$ADM2_PT[which(area_municipality$harv_area_frac > 1, arr.ind=TRUE)]) +#summary(area_municipality$harv_area_frac, na.rm = TRUE ) +#hist(area_municipality$harv_area_frac) + +# Check summary data for harvest area +# summary(shp.soy.has$Valor ) +# hist(area_municipality$Valor) +# test_2 <-area_municipality$Valor +# test_3 <- with(area_municipality,replace(Valor, harv_area_frac < 0.01, NA)) +# summary(test_2 ) +# summary(test_3 ) +# hist(test_2) +# hist(test_3) +# all.equal(test_2,test_3) +# sum(test_2 > 0, na.rm=T) - sum(test_3 > 0, na.rm=T) # number of cases that disappeared + +shp.soy.has$harv_area_frac <- area_municipality$harv_area_frac +shp.soy.has$Valor <- with(area_municipality,replace(Valor, harv_area_frac < 0.01, NA)) +sum(shp.soy.has$Valor[shp.soy.has$harv_area_frac < 0.01], na.rm=T) # IT should be zero to show there is no value below 1% + +summary(shp.soy.has$Valor ) +summary(shp.soy.has$harv_area_frac ) + + +shp.soy.has.2009 <- shp.soy.has %>% filter(Ano == 2016) +# plot(shp.soy.has.2009['Valor']) + +# Save shapefile with the cells only with > 1% + +#### +extent_br = c(-180, 180, -90, 90) +resolution_br = 0.5 +ncol_br = (extent_br[2]-extent_br[1])/resolution_br +nrows_br = (extent_br[4]-extent_br[3])/resolution_br + +# Raster creation +baserast <- rast(nrows=nrows_br, ncol=ncol_br, + extent= extent_br, + crs="+proj=longlat +datum=WGS84") + +# rasters <- rast(lapply(1980:2016, +# function(x) +# rasterize(vect(shp.soy.has %>% +# filter(Ano==x)),baserast,"Valor", fun = 'mean'))) +# names(rasters) <- 1980:2016 +# varnames(rasters) <- paste0("soy_yield_",1980:2016) +# +# writeRaster(rasters,"soy_harvest_area_city_1980_2016_1prc_05x05.tiff",overwrite=TRUE) +# +# plot(rasters) + + + +###### Density for harvest area +r <- baserast +shp.soy.yld.br_subset <- subset(shp.soy.has, Ano < 2017 ) +v <- vect(shp.soy.yld.br_subset) +ra <- cellSize(r, unit="ha") +e <- expanse(v, unit="ha") +v$density <- v$Valor / e + +years <- str_sort(unique(v$Ano)) +out <- list() +for (i in 1:length(years)) { + vv <- v[v$Ano == years[i], ] + x <- rasterize(vv, r, "density") + out[[i]] <- x * ra +} +out <- rast(out) +names(out) <- years + +writeRaster(out,"soy_harvest_area_br_1980_2016_05x05_density_03.tif", overwrite=TRUE) + + diff --git a/code/create_raster_soy_yield_br.R b/code/create_raster_soy_yield_br.R new file mode 100644 index 0000000..319d4a9 --- /dev/null +++ b/code/create_raster_soy_yield_br.R @@ -0,0 +1,195 @@ +library(sidrar) +library(terra) +library(tidyr) +library(dplyr) +library(sf) +library(data.table) +setwd("~/PhD/paper_hybrid_agri/data") +dir_crop_ibge <- ("C:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data") + +init_date <- 1975 +end_date <- 2016 + +# period.list <- list(1975:1979) +# # info_sidra(1612)[[3]] +# var_list <- c(112) +# +# mun.soy.yld <- lapply(period.list, function(x) +# get_sidra(1612, variable = var_list, +# period = as.character(x), +# geo = "City", +# classific = c("c81"), +# category = list(2713))) +# +# test_2 <- do.call(rbind,mun.soy.yld) +# write.csv(test_2,"soy_yield_munic_br_1975_1978.csv") + +mun.soy.yld1 <- read.csv("soy_yield_munic_br_1975_1978.csv") +mun.soy.yld1 <- subset(mun.soy.yld1, Ano..Código. < 1979) + +mun.soy.yld2 <- read.csv("soy_yield_munic_br.csv") +mun.soy.yld <- rbind(mun.soy.yld1, mun.soy.yld2) + +# mun.soy.yld3 <- list.files(dir_crop_ibge, "soy_yield_munic_br", full.names = TRUE) %>% +# map_dfr(read.csv) + +mun.soy.yld$Valor <- mun.soy.yld$Valor / 1000 + +mun.soy.yld <- mun.soy.yld %>% rename(ADM2_PCODE = 'Município..Código.') %>% + dplyr::select(Ano,ADM2_PCODE,Valor) + +mun.soy.yld <- dplyr::arrange(mun.soy.yld , ADM2_PCODE, Ano) + + +BRmun <- st_read("GIS/bra_admbnda_adm2_ibge_2020.shp") %>% + mutate(ADM2_PCODE = substr(ADM2_PCODE,3,9)) %>% + mutate_at(c('ADM2_PCODE'), as.numeric) + +st_crs(BRmun) <- 4326 +BRmun <- st_transform(BRmun, crs = 4326) + +shp.soy.yld.br <- left_join(BRmun,mun.soy.yld) %>% drop_na() + +# RASTER BASE EXTENTION +extent_br = c(-180, 180, -90, 90) +resolution_br = 0.05 +ncol_br = (extent_br[2]-extent_br[1])/resolution_br +nrows_br = (extent_br[4]-extent_br[3])/resolution_br + +# Raster creation +baserast <- rast(nrows=nrows_br, ncol=ncol_br, + extent= extent_br, + crs="+proj=longlat +datum=WGS84", + vals=NA) + +rasters <- rast(lapply(init_date:end_date, + function(x) + rasterize(vect(shp.soy.yld.br %>% + filter(Ano==x)),baserast,"Valor"))) +names(rasters) <- init_date:end_date +varnames(rasters) <- paste0("soy_yield_",init_date:end_date) + + +writeRaster(rasters,"soy_yield_1975_2016_005.tif", overwrite = TRUE) + + + +### Add the conditions for harvested area > 1% ################################# +test <- shp.soy.yld.br %>% st_drop_geometry() +test_2 <- shp.soy.has %>% st_drop_geometry() +remove_cities <- anti_join(test_2, test, by = c('ADM2_PCODE','Ano','ADM1_PCODE')) + +shp.soy.yld3 <- shp.soy.yld.br + +hist(shp.soy.yld3$Valor) +summary(shp.soy.yld3$Valor) + +shp.soy.yld3$harv_area_frac <- area_municipality$harv_area_frac +sum(shp.soy.yld3$harv_area_frac < 0.01, na.rm=T) +sum(shp.soy.yld3$Valor[shp.soy.yld3$harv_area_frac < 0.01], na.rm=T) + +shp.soy.yld3$Valor[shp.soy.yld3$harv_area_frac < 0.01] <- NA + +sum(shp.soy.yld3$harv_area_frac < 0.01, na.rm=T) +sum(shp.soy.yld3$Valor[shp.soy.yld3$harv_area_frac < 0.01], na.rm=T) + +extent_br = c(-180, 180, -90, 90) +resolution_br = 0.05 +ncol_br = (extent_br[2]-extent_br[1])/resolution_br +nrows_br = (extent_br[4]-extent_br[3])/resolution_br + +# Raster creation +baserast <- rast(nrows=nrows_br, ncol=ncol_br, + extent= extent_br, + crs="+proj=longlat +datum=WGS84", + vals=NA) + +rasters <- rast(lapply(init_date:end_date, + function(x) + rasterize(vect(shp.soy.yld3 %>% + filter(Ano==x)),baserast,"Valor"))) +names(rasters) <- init_date:end_date +varnames(rasters) <- paste0("soy_yield_",init_date:end_date) + +writeRaster(rasters,"soy_yield_1975_2016_005_1prc.tif", overwrite=TRUE) + + +### Filters +# (1) counties with more than 3 years (???4) of repeated values were considered to contain artificial data and were excluded from the analysis; +# (2) only counties with more than 4 years (???5) of consecutive records were selected for the analysis, + +### Remove counties that have 6 or more times repeated the same value ################################# +df_repeat <- shp.soy.yld3 %>% + count(Valor, ADM2_PCODE) %>% + filter(!is.na(Valor)) %>% + group_by(ADM2_PCODE) %>% + summarise(Value = max(n, na.rm = TRUE)) %>% + mutate(repeat_mask = ifelse(Value >= 6, TRUE, FALSE)) + +### Remove counties where there are 4 or more consecutive identical values ################################# +df_repeat2 <- shp.soy.yld3 %>% + filter(!is.na(Valor)) %>% + group_by(ADM2_PCODE) %>% + summarise(max_consecutive_values = max(rle(Valor)$lengths, na.rm = TRUE)) %>% + mutate(repeat_mask = ifelse(max_consecutive_values >= 4, TRUE, FALSE)) + +### Remove counties with 10 or less years ################################# +df_minimum_years <- shp.soy.yld3 %>% + filter(!is.na(Valor)) %>% + group_by(ADM2_PCODE) %>% + count(duration_year = max(as.integer(Ano))-min(as.integer(Ano))) %>% + mutate(time_mask = ifelse(n >= 10, FALSE, TRUE)) +hist(df_minimum_years$n) + +# Check amount of counties being removed for each criterium +sum(df_repeat$repeat_mask == TRUE) +sum(df_repeat2$repeat_mask == TRUE) +sum(df_minimum_years$time_mask == TRUE) + +df_repeat <- data.frame(df_repeat$ADM2_PCODE,df_repeat$repeat_mask) +names(df_repeat) <- c("ADM2_PCODE", "repeat_mask") + +df_repeat2 <- data.frame(df_repeat2$ADM2_PCODE,df_repeat2$repeat_mask) +names(df_repeat2) <- c("ADM2_PCODE", "consecutive_mask") + +df_minimum_years <- data.frame(df_minimum_years$ADM2_PCODE,df_minimum_years$time_mask) +names(df_minimum_years) <- c("ADM2_PCODE", "minimum_time_mask") + +df_all_filters <- data.frame(df_repeat$ADM2_PCODE, df_repeat$repeat_mask, df_repeat2$minimum_time_mask,df_minimum_years$minimum_time_mask) +names(df_all_filters) <- c("ADM2_PCODE","repeat_mask","consecutive_mask", "minimum_time_mask") + +# df_soy <- data.frame(shp.soy.yld3$ADM2_PCODE,shp.soy.yld3$Valor) +# names(df_soy) <- c("ADM2_PCODE", "Valor") +# +# df_repeat3 <- df_soy %>% +# filter(!is.na(Valor)) %>% +# group_by(ADM2_PCODE) %>% +# mutate(grp = rleid(Valor)) %>% +# count(grp) %>% +# summarise(max_consecutive_values = max(n)) + +shp.soy.yld4 <- left_join(shp.soy.yld3,df_all_filters) +shp.soy.yld4$Valor[shp.soy.yld4$repeat_mask == TRUE ] <- NA +shp.soy.yld4$Valor[shp.soy.yld4$consecutive_mask == TRUE ] <- NA +shp.soy.yld4$Valor[shp.soy.yld4$minimum_time_mask == TRUE ] <- NA +hist(shp.soy.yld3$Valor) +hist(shp.soy.yld4$Valor) +summary(shp.soy.yld3$Valor) +summary(shp.soy.yld4$Valor) + + +### Create Raster +baserast <- rast(nrows=4860, ncol=5040, + extent= c(-74.25, -32.25, -34.25, 6.25), + crs="+proj=longlat +datum=WGS84", + vals=NA) + +rasters <- rast(lapply(1980:2016, + function(x) + rasterize(vect(shp.soy.yld4 %>% + filter(Ano==x)),baserast,"Valor"))) +names(rasters) <- 1980:2016 +varnames(rasters) <- paste0("soy_yield_",1980:2016) + +writeRaster(rasters,"soy_yield_1980_2016_all_filters.tiff") + diff --git a/code/create_raster_soy_yield_us.R b/code/create_raster_soy_yield_us.R new file mode 100644 index 0000000..9faf032 --- /dev/null +++ b/code/create_raster_soy_yield_us.R @@ -0,0 +1,129 @@ +library(terra) +library(tidyr) +library(dplyr) +library(tidyverse) +library(USAboundaries) +library(stars) + +# FOR HARVEST YIELD_us +setwd("~/PhD/paper_hybrid_agri/data") + +crop_dir <- ("C:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/Paper_drought/data/usda_data") + +# usda_yields <- read.csv("soy_usda_obs_1960_2019.csv") %>% dplyr::select(County,State,Year,Value) + +usda_yields <- list.files(crop_dir, "usda", full.names = TRUE) %>% + map_dfr(read.csv) %>% + arrange(Year) %>% dplyr::select(County,State,Year,Value) + +usda_yields$Value <- usda_yields$Value * 0.06725106937 + +usda_yields <- subset(usda_yields, (County !="SEVIER" | State != 'ARKANSAS' | Year != '1989')) + +us_spatial_sf<- function(x) { + #x is either us_counties or us_states from usaboundaries package + x %>% + filter(!state_name %in% + c("Alaska","Puerto Rico","Hawaii","Washington", + "Oregon","California","Nevada","Idaho","Utah", + "Arizona","New Mexico","Colorado","Wyoming","Montana"))%>% + rename(County = name, + State = state_name) %>% + dplyr::select(County,State,aland,geometry) %>% + mutate_at(vars(State,County), ~toupper(.)) +} + +us_counties_subset <- us_counties() %>% dplyr::select(c(6,9,11,15)) + +usda_spatial_full <- us_counties_subset %>% + us_spatial_sf %>% + left_join(usda_yields, by = c("County", "State")) %>% + drop_na() + +usda_spatial_yield <-usda_spatial_full %>% + dplyr::select(County,State,Year,Value,geometry) + +usda_spatial_2010 <- usda_spatial_yield[usda_spatial_yield$Year == 1980,] +usda_spatial_2010["Value"] %>% plot() + +extent_us = c(-180, 180, -90, 90) +resolution_us = 0.05 +ncol_us = (extent_us[2]-extent_us[1])/resolution_us +nrows_us = (extent_us[4]-extent_us[3])/resolution_us + +# Raster creation +baserast <- rast(nrows=nrows_us, ncol=ncol_us, + extent= extent_us, # CHANGE TO US + crs="+proj=longlat +datum=WGS84", + vals=NA) + +rasters <- rast(lapply(1975:2020, + function(x) + rasterize(vect(usda_spatial_yield %>% + filter(Year==x)),baserast,"Value"))) +names(rasters) <- 1975:2020 +varnames(rasters) <- paste0("soy_yield_usda",1975:2020) + +writeRaster(rasters,"soy_yields_US_all_1975_2020_001.tif",overwrite=TRUE) + +######################################################################################### +#### ADD filter for 1% of removal - needs file create_raster_harvest_area_us.R loaded +# setwd("~/PhD/paper_hybrid_agri/data") + +# Calculate harvest area % +total_area <- usda_spatial_full %>% + dplyr::select(County,State,Year,aland,geometry) %>% dplyr::mutate(aland = aland/10000) + + +usda_spatial_yield$ratio_area <- usda_spatial_area$Value / total_area$aland + +sum(usda_spatial_yield$ratio_area < 0.01, na.rm=T) +sum(usda_spatial_yield$Value[usda_spatial_yield$ratio_area < 0.01], na.rm=T) + +usda_spatial_yield$Value <- with(usda_spatial_yield,replace(Value, ratio_area < 0.01, NA)) +sum(usda_spatial_yield$Value[usda_spatial_yield$ratio_area < 0.01], na.rm=T) + + +usda_spatial_2010 <- usda_spatial_yield[usda_spatial_yield$Year == 1980,] +usda_spatial_2010["Value"] %>% plot() + +rasters <- rast(lapply(1975:2020, + function(x) + rasterize(vect(usda_spatial_yield %>% + filter(Year==x)),baserast,"Value"))) +names(rasters) <- 1975:2020 +varnames(rasters) <- paste0("soy_yield_usda",1975:2020) + +writeRaster(rasters,"soy_yields_US_all_1975_2020_1prc_002.tif",overwrite=TRUE) + + + + +# +# ################################################################################ +# # VECTORISE EPIC OUTPUT +# rdir <- ("C:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/Paper_drought/data") +# epic_yield_us = read_ncdf(paste(rdir,"ACY_gswp3-w5e5_obsclim_2015soc_default_soy_noirr.nc",sep="")) +# +# +# us_climate <- us_climate_xy%>% +# st_as_sf(x = ., coords = c("x", "y"), +# crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") +# +# +# #====================================================================== +# #Aggregate grid based climate information to county scale +# df_point_list <- split(dplyr::select(us_climate, -Year), +# us_climate$Year) +# #Split sf usda_resid country by Years +# df_poly_list <- split(usda_resid_county, usda_resid_county$Year) +# #join yield and climate sf, na are reproduced as climate_df is filtered +# full_Sf_yield_climate <- map2_dfr(df_poly_list, df_point_list, +# ~ .x %>% +# st_join(.y, left =FALSE) %>% +# group_by(County,Year,State) %>% +# mutate_at(vars(-geometry,-State,-County,-Year),~(mean(.,na.rm = TRUE)))) +# #do you want a simple spatial averages? +# #summarize_at(vars(-geometry),~(mean(.,na.rm = TRUE)))) +# #### +# diff --git a/code/ga_sim_plots.py b/code/ga_sim_plots.py index 093ff37..2225d87 100644 --- a/code/ga_sim_plots.py +++ b/code/ga_sim_plots.py @@ -13,7 +13,9 @@ """ import os -os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/Paper_drought/data') +os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data/ga_simulations') + +# os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/Paper_drought/data') import xarray as xr import numpy as np @@ -25,11 +27,7 @@ import cartopy.io.shapereader as shpreader from scipy import signal import seaborn as sns -from mask_shape_border import mask_shape_border -from failure_probability import feature_importance_selection, failure_probability -from stochastic_optimization_Algorithm import stochastic_optimization_Algorithm -from shap_prop import shap_prop -from bias_correction_masked import * +# from mask_shape_border import mask_shape_border import matplotlib as mpl import pickle @@ -177,6 +175,8 @@ def shock_conversion(prod_soya_br_subset, chosen_model = None, co2_scenario = No plt.legend( title="",loc=3, fontsize='small', fancybox=True) plt.show() +df_bar_br_us.to_csv("aggregated_shocks_soybean_br_us_gw.csv") + #%% df_soy_glob_subset_xprp = df_soy_glob[ (df_soy_glob['VAR_ID'] == 'XPRP') & (df_soy_glob['VAR_UNIT'] == 'fm t/ha') & (df_soy_glob['MACROSCEN'] == 'SSP2') ] @@ -269,7 +269,6 @@ def shock_conversion(prod_soya_br_subset, chosen_model = None, co2_scenario = No def plot_shocks_socio(Item_sel, Var_id_sel, Region = 'EU'): df_europe_imp_scen = df_europe_imp_subset[ (df_europe_imp_subset['VAR_ID'] == Var_id_sel) & (df_europe_imp_subset['Item'] == Item_sel) & (df_europe_imp_subset['ANYREGION'] == Region)] df_europe_imp_ref_sub = df_europe_imp_ref[(df_europe_imp_ref['VAR_ID'] == Var_id_sel) & (df_europe_imp_ref['Item'] == Item_sel)& (df_europe_imp_ref['ANYREGION'] == Region)] - print(df_europe_imp_scen) # BAR PLOTS EUROPE plt.figure(figsize = (7,5),dpi=300) sns.barplot(x=df_europe_imp_scen["GW degree"], y=df_europe_imp_scen["VALUE"], hue = df_europe_imp_scen['BIOENSCEN'], order = ['2.5°C', '3°C']) @@ -278,6 +277,7 @@ def plot_shocks_socio(Item_sel, Var_id_sel, Region = 'EU'): plt.ylabel('Shock (%)') plt.legend() plt.show() + df_europe_imp_scen.to_csv(f'output_globiom_gw_{Item_sel}_{Var_id_sel}_{Region}.csv') plot_shocks_socio('Soya','CONS') plot_shocks_socio('Soya','XPRP') @@ -304,7 +304,7 @@ def plot_shocks_socio(Item_sel, Var_id_sel, Region = 'EU'): #%% Change visualization to degree of GW import glob -def open_csv_timeseries(path= "projections_global_mean/ipsl-cm6a-lr_r1i1p1f1_w5e5_ssp126_tas_*.csv", pre_ind_tmep = 13.8): +def open_csv_timeseries(path, pre_ind_tmep = 13.8): files = glob.glob(path) df = [] for f in files: @@ -317,16 +317,16 @@ def open_csv_timeseries(path= "projections_global_mean/ipsl-cm6a-lr_r1i1p1f1_w5e return df_2 -df_ipsl_26 = open_csv_timeseries("../projections_global_mean/ipsl-cm6a-lr_r1i1p1f1_w5e5_ssp126_tas_*.csv") -df_ipsl_85 = open_csv_timeseries("../projections_global_mean/ipsl-cm6a-lr_r1i1p1f1_w5e5_ssp585_tas_*.csv") +df_ipsl_26 = open_csv_timeseries("projections_global_mean/ipsl-cm6a-lr_r1i1p1f1_w5e5_ssp126_tas_*.csv") +df_ipsl_85 = open_csv_timeseries("projections_global_mean/ipsl-cm6a-lr_r1i1p1f1_w5e5_ssp585_tas_*.csv") -df_ukesm_26 = open_csv_timeseries("../projections_global_mean/ukesm1-0-ll_r1i1p1f2_w5e5_ssp126_tas_*.csv") -df_ukesm_85 = open_csv_timeseries("../projections_global_mean/ukesm1-0-ll_r1i1p1f2_w5e5_ssp585_tas_*.csv") +df_ukesm_26 = open_csv_timeseries("projections_global_mean/ukesm1-0-ll_r1i1p1f2_w5e5_ssp126_tas_*.csv") +df_ukesm_85 = open_csv_timeseries("projections_global_mean/ukesm1-0-ll_r1i1p1f2_w5e5_ssp585_tas_*.csv") -df_gfdl_26 = open_csv_timeseries("../projections_global_mean/gfdl-esm4_r1i1p1f1_w5e5_ssp126_tas_*.csv") -df_gfdl_85 = open_csv_timeseries("../projections_global_mean/gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_*.csv") +df_gfdl_26 = open_csv_timeseries("projections_global_mean/gfdl-esm4_r1i1p1f1_w5e5_ssp126_tas_*.csv") +df_gfdl_85 = open_csv_timeseries("projections_global_mean/gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_*.csv") df_ipsl_26.groupby(df_ipsl_26.index.year)['tas'].transform('mean').plot() df_ipsl_85.groupby(df_ipsl_85.index.year)['tas'].transform('mean').plot() diff --git a/code/globiom_grid_fun_am.R b/code/globiom_grid_fun_am.R new file mode 100644 index 0000000..6cf099f --- /dev/null +++ b/code/globiom_grid_fun_am.R @@ -0,0 +1,323 @@ +library(ncdf4) # package for netcdf manipulation +library(raster) # package for raster manipulation +library(rgdal) # package for geospatial analysis +library(ggplot2) # package for plotting +library(terra) +library(tidyr) +library(dplyr) +library(sf) +library(sp) + +rdir <- ("C:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data/output_shocks_am/scenarios_forum") +setwd(rdir) + +globiom_grid_am <- vect("../../soy_americas_harvest_area.shp") + + +### MAIN CODE +create_shifters_globiom <- function(model, co2_scenario, ssp, years, three_models = FALSE) { + + # Load .nc file + nc_file <- raster(paste('hybrid_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.nc', sep = "")) + + writeRaster(x = nc_file, filename = paste('hybrid_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.tif', sep = ""), format = 'GTiff', overwrite = TRUE) + + raster_2012_hybrid <- rast(paste('hybrid_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.tif', sep = "")) + names(raster_2012_hybrid) <- 'hybrid_shift' + + # If only using Hybrid + if (three_models == FALSE) { + v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) + + output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) + + globiom_grid_am$shift_hybrid <- output_hybrid + + shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) + names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') + shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + } + + + if (three_models == TRUE) { + writeRaster(x =raster(paste('epic_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.nc', sep = "")), filename = paste('epic_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.tif', sep = ""), format = 'GTiff', overwrite = TRUE) + + writeRaster(x =raster(paste('clim_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.nc', sep = "")), filename = paste('clim_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.tif', sep = ""), format = 'GTiff', overwrite = TRUE) + + raster_2012_epic <- rast(paste('epic_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.tif', sep = "")) + names(raster_2012_epic) <- 'epic_shift' + + raster_2012_clim <- rast(paste('clim_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.tif', sep = "")) + names(raster_2012_clim) <- 'clim_shift' + + all_rasters <- c(raster_2012_hybrid,raster_2012_epic,raster_2012_clim ) + + v_all <- terra::extract(all_rasters, globiom_grid_am, fun = mean) + + output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) + output_epic = unlist(lapply(v_all$epic_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) + output_clim = unlist(lapply(v_all$clim_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) + + globiom_grid_am$shift_hybrid <- output_hybrid + globiom_grid_am$shift_epic <- output_epic + globiom_grid_am$shift_clim <- output_clim + + shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid, globiom_grid_am$shift_epic, globiom_grid_am$shift_clim) + names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid', 'epic', 'clim') + shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + } + + # Save shifters to csv + write.csv(shift_2012_soybean_am_sub, paste('output_csv/yield_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_',years,'.csv', sep = ""), row.names = FALSE) + +} + +##### MERGE CSV FILES +merge_shifters <- function(model,co2_scenario) { + + group_per_time <- function(model, co2_scenario, ssp) { + + csv_shift_2017 <- read.csv(paste('output_csv/yield_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_2017-2044.csv', sep = "")) + csv_shift_2044 <- read.csv(paste('output_csv/yield_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_2044-2071.csv', sep = "")) + csv_shift_2071 <- read.csv(paste('output_csv/yield_', model, '_',ssp,'_',co2_scenario,'_yield_soybean_shift_2071-2098.csv', sep = "")) + + csv_shift_2017$ALLYEAR = 2030 + csv_shift_2044$ALLYEAR = 2050 + csv_shift_2071$ALLYEAR = 2070 + + csv_combined <- bind_rows(csv_shift_2017,csv_shift_2044,csv_shift_2071 ) + + csv_combined <- csv_combined %>% select(ALLYEAR, everything()) + + if(co2_scenario == '2015co2' && ssp == 'ssp585' ){ + csv_combined$rcp = 'noC8p5' + } else if (co2_scenario == 'default' && ssp == 'ssp585') { + csv_combined$rcp = 'rcp8p5' + } else if (co2_scenario == '2015co2' && ssp == 'ssp126') { + csv_combined$rcp = 'noC2p6' + } else if (co2_scenario == 'default' && ssp == 'ssp126') { + csv_combined$rcp = 'rcp2p6' + } + + write.csv(csv_combined, paste('output_csv/yield_', model, '_',ssp,'_',co2_scenario,'_soybean_shift_2015-2100.csv', sep = ""), row.names = FALSE) + + } + + group_per_time(model, co2_scenario, 'ssp585') + group_per_time(model, co2_scenario, 'ssp126') + + ###### COMBINE TWO RCPs + csv_shift_85 <- read.csv(paste('output_csv/yield_', model, '_ssp585_',co2_scenario,'_soybean_shift_2015-2100.csv', sep = "")) + csv_shift_26 <- read.csv(paste('output_csv/yield_', model, '_ssp126_',co2_scenario,'_soybean_shift_2015-2100.csv', sep = "")) + + csv_comb_rcps <- bind_rows(csv_shift_26,csv_shift_85 ) + csv_comb_rcps <- csv_comb_rcps %>% select(rcp, everything()) + + indir <- paste(rdir,"/results_shifters_am_",model,'_ssp126_585',"/",sep="") + dir.create(indir) + + write.csv(csv_comb_rcps, paste(indir, 'yield_', model, '_ssp126_585_',co2_scenario,'_soybean_shift_2015-2100.csv', sep = ""), row.names = FALSE) + +} + +merge_models <- function(model1,model2,model3, co2_scenario) { + csv_shift_gfdl <- read.csv(paste(rdir,"/results_shifters_am_",model1,'_ssp126_585',"/",'yield_', model1, '_ssp126_585_',co2_scenario,'_soybean_shift_2015-2100.csv',sep="")) + csv_shift_ukesm <- read.csv(paste(rdir,"/results_shifters_am_",model2,'_ssp126_585',"/",'yield_', model2, '_ssp126_585_',co2_scenario,'_soybean_shift_2015-2100.csv',sep="")) + csv_shift_ipsl <- read.csv(paste(rdir,"/results_shifters_am_",model3,'_ssp126_585',"/",'yield_', model3, '_ssp126_585_',co2_scenario,'_soybean_shift_2015-2100.csv',sep="")) + + csv_shift_gfdl$GCM = 'gfdl-esm4' + csv_shift_ukesm$GCM = 'ukesm1-0-ll' + csv_shift_ipsl$GCM = 'ipsl-cm6a-lr' + + csv_combined <- bind_rows(csv_shift_gfdl,csv_shift_ukesm,csv_shift_ipsl ) + csv_combined <- csv_combined %>% select(GCM, everything()) + + indir <- paste(rdir,"/results_shifters_am_all_gcms_rcps","/",sep="") + dir.create(indir) + + write.csv(csv_combined, paste(indir, 'yield_all_models_rcps_soybean_shift_2015-2100.csv', sep = ""), row.names = FALSE) +} + +# +# # 2015co2 +# create_shifters_globiom('gfdl-esm4', '2015co2', 'ssp126' ,'2017-2044') +# create_shifters_globiom('gfdl-esm4', '2015co2','ssp126' ,'2044-2071') +# create_shifters_globiom('gfdl-esm4', '2015co2','ssp126' ,'2071-2098') +# +# create_shifters_globiom('gfdl-esm4', '2015co2','ssp585' ,'2017-2044') +# create_shifters_globiom('gfdl-esm4', '2015co2','ssp585' ,'2044-2071') +# create_shifters_globiom('gfdl-esm4', '2015co2','ssp585' ,'2071-2098') +# +# +# create_shifters_globiom('ukesm1-0-ll', '2015co2', 'ssp126' ,'2017-2044') +# create_shifters_globiom('ukesm1-0-ll', '2015co2','ssp126' ,'2044-2071') +# create_shifters_globiom('ukesm1-0-ll', '2015co2','ssp126' ,'2071-2098') +# +# create_shifters_globiom('ukesm1-0-ll', '2015co2','ssp585' ,'2017-2044') +# create_shifters_globiom('ukesm1-0-ll', '2015co2','ssp585' ,'2044-2071') +# create_shifters_globiom('ukesm1-0-ll', '2015co2','ssp585' ,'2071-2098') +# +# +# create_shifters_globiom('ipsl-cm6a-lr', '2015co2','ssp126' ,'2017-2044') +# create_shifters_globiom('ipsl-cm6a-lr', '2015co2','ssp126' ,'2044-2071') +# create_shifters_globiom('ipsl-cm6a-lr', '2015co2','ssp126' ,'2071-2098') +# +# create_shifters_globiom('ipsl-cm6a-lr', '2015co2','ssp585' ,'2017-2044') +# create_shifters_globiom('ipsl-cm6a-lr', '2015co2','ssp585' ,'2044-2071') +# create_shifters_globiom('ipsl-cm6a-lr', '2015co2','ssp585' ,'2071-2098') + + +# DEFAULT co2 +create_shifters_globiom('gfdl-esm4', 'default','ssp126' ,'2017-2044') +create_shifters_globiom('gfdl-esm4', 'default','ssp126' ,'2044-2071') +create_shifters_globiom('gfdl-esm4', 'default','ssp126' ,'2071-2098') + +create_shifters_globiom('gfdl-esm4', 'default','ssp585' ,'2017-2044') +create_shifters_globiom('gfdl-esm4', 'default','ssp585' ,'2044-2071') +create_shifters_globiom('gfdl-esm4', 'default','ssp585' ,'2071-2098') + +create_shifters_globiom('ukesm1-0-ll', 'default','ssp126' ,'2017-2044') +create_shifters_globiom('ukesm1-0-ll', 'default','ssp126' ,'2044-2071') +create_shifters_globiom('ukesm1-0-ll', 'default','ssp126' ,'2071-2098') + +create_shifters_globiom('ukesm1-0-ll', 'default','ssp585' ,'2017-2044') +create_shifters_globiom('ukesm1-0-ll', 'default','ssp585' ,'2044-2071') +create_shifters_globiom('ukesm1-0-ll', 'default','ssp585' ,'2071-2098') + +create_shifters_globiom('ipsl-cm6a-lr', 'default','ssp126' ,'2017-2044') +create_shifters_globiom('ipsl-cm6a-lr', 'default','ssp126' ,'2044-2071') +create_shifters_globiom('ipsl-cm6a-lr', 'default','ssp126' ,'2071-2098') + +create_shifters_globiom('ipsl-cm6a-lr', 'default','ssp585' ,'2017-2044') +create_shifters_globiom('ipsl-cm6a-lr', 'default','ssp585' ,'2044-2071') +create_shifters_globiom('ipsl-cm6a-lr', 'default','ssp585' ,'2071-2098') + +# +# ### FINAL SHIFTERS FOR MODEL AND CO2 Scenarios +# merge_shifters('gfdl-esm4', '2015co2') +# merge_shifters('ukesm1-0-ll', '2015co2') +# merge_shifters('ipsl-cm6a-lr', '2015co2') + +merge_shifters('gfdl-esm4', 'default') +merge_shifters('ukesm1-0-ll', 'default') +merge_shifters('ipsl-cm6a-lr', 'default') + + +# MERGE MODEL +merge_models('gfdl-esm4','ukesm1-0-ll','ipsl-cm6a-lr', 'default') + + +############### Historical run - 2012 +# load file and concvert to tiff. +nc_file_hist <- raster(paste('../../output_models_am/shifters_2012_hybrid_climatology_syn_2050area.nc', sep = "")) +writeRaster(x = nc_file_hist, filename = paste('../../output_models_am/shifters_2012_hybrid_climatology_syn_2050area','.tif', sep = ""), format = 'GTiff', overwrite = TRUE) +raster_2012_hybrid <- rast(paste('../../output_models_am/shifters_2012_hybrid_climatology_syn_2050area','.tif', sep = "")) +names(raster_2012_hybrid) <- 'hybrid_shift' +# Process the data and convert to csv according to globiom grid +v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) +output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) +globiom_grid_am$shift_hybrid <- output_hybrid +shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) +names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') +shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + +# Save shifters to csv +write.csv(shift_2012_soybean_am_sub, paste('output_csv/shifters_2012_hybrid_climatology_syn_2050area','.csv', sep = ""), row.names = FALSE) + + +## ESTHER - future shifters +nc_file_hist <- raster(paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp126_default_yield_soybean_shift_2044-2071.nc', sep = "")) +writeRaster(x = nc_file_hist, filename = paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp126_default_yield_soybean_shift_2044-2071','.tif', sep = ""), format = 'GTiff', overwrite = TRUE) +raster_2012_hybrid <- rast(paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp126_default_yield_soybean_shift_2044-2071','.tif', sep = "")) +names(raster_2012_hybrid) <- 'hybrid_shift' +# Process the data and convert to csv according to globiom grid +v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) +output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) +globiom_grid_am$shift_hybrid <- output_hybrid +shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) +names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') +shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) +# Save shifters to csv +write.csv(shift_2012_soybean_am_sub, paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp126_default_yield_soybean_shift_2044-2071','.csv', sep = ""), row.names = FALSE) + +nc_file_hist <- raster(paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp585_default_yield_soybean_shift_2044-2071.nc', sep = "")) +writeRaster(x = nc_file_hist, filename = paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp585_default_yield_soybean_shift_2044-2071','.tif', sep = ""), format = 'GTiff', overwrite = TRUE) +raster_2012_hybrid <- rast(paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp585_default_yield_soybean_shift_2044-2071','.tif', sep = "")) +names(raster_2012_hybrid) <- 'hybrid_shift' +# Process the data and convert to csv according to globiom grid +v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) +output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) +globiom_grid_am$shift_hybrid <- output_hybrid +shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) +names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') +shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) +# Save shifters to csv +write.csv(shift_2012_soybean_am_sub, paste('../../output_shocks_am/esther_paper/hybrid_ukesm1-0-ll_ssp585_default_yield_soybean_shift_2044-2071','.csv', sep = ""), row.names = FALSE) + + + + + +############### SAVE GLOBIOM GRID +globiom_grid_am <- vect("../../soy_americas_harvest_area.shp") +# load file and concvert to tiff. +nc_file_hist <- raster(paste('../../soy_harvest_spam_native_05x05.nc', sep = "")) +writeRaster(x = nc_file_hist, filename = paste('../../soy_harvest_spam_native_05x05','.tif', sep = ""), format = 'GTiff', overwrite = TRUE) + +raster_2012_hybrid <- rast(paste('../../soy_harvest_spam_native','.tif', sep = "")) +names(raster_2012_hybrid) <- 'hybrid_shift' + +# Process the data and convert to csv according to globiom grid +v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = sum) +output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) +globiom_grid_am$shift_hybrid <- output_hybrid +shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) +names(shift_2012_soybean_am) <- c('ALLCOLROW','harvest_area') +shift_2012_soybean_am_sub3 <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + +# Save shifters to csv +write.csv(shift_2012_soybean_am_sub, paste('../../soy_harvest_spam_globiom','.csv', sep = ""), row.names = FALSE) + + + +# Sendin the minimum shocks in the entire series to Esther +create_shifters_globiom_min <- function(input) { + nc_file_hist <- raster(paste(input,'.nc', sep = "")) + writeRaster(x = nc_file_hist, filename = paste(input,'.tif', sep = ""), format = 'GTiff', overwrite = TRUE) + raster_2012_hybrid <- rast(paste(input,'.tif', sep = "")) + names(raster_2012_hybrid) <- 'hybrid_shift' + # Process the data and convert to csv according to globiom grid + v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) + output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) + globiom_grid_am$shift_hybrid <- output_hybrid + shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) + names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') + shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + + str_parts = strsplit(input, "/") + # Save shifters to csv + write.csv(shift_2012_soybean_am_sub, paste(str_parts[[1]][1], '/', str_parts[[1]][2], '/output_csv/', str_parts[[1]][3],'.csv', sep = ""), row.names = FALSE)} + + +create_shifters_globiom_min('../esther_paper/hybrid_gfdl-esm4_ssp126_default_yield_soybean_min_shift_2015_2098') +create_shifters_globiom_min('../esther_paper/hybrid_gfdl-esm4_ssp585_default_yield_soybean_min_shift_2015_2098') +create_shifters_globiom_min('../esther_paper/hybrid_ipsl-cm6a-lr_ssp126_default_yield_soybean_min_shift_2015_2098') +create_shifters_globiom_min('../esther_paper/hybrid_ipsl-cm6a-lr_ssp585_default_yield_soybean_min_shift_2015_2098') +create_shifters_globiom_min('../esther_paper/hybrid_ukesm1-0-ll_ssp126_default_yield_soybean_min_shift_2015_2098') +create_shifters_globiom_min('../esther_paper/hybrid_ukesm1-0-ll_ssp585_default_yield_soybean_min_shift_2015_2098') + + +nc_file_hist <- raster(paste('../esther_paper/hybrid_gfdl-esm4_ssp126_default_yield_soybean_min_shift_2015_2098', '.nc', sep = "")) +writeRaster(x = nc_file_hist, filename = paste('../esther_paper/hybrid_gfdl-esm4_ssp126_default_yield_soybean_min_shift_2015_2098','.tif', sep = ""), format = 'GTiff', overwrite = TRUE) +raster_2012_hybrid <- rast(paste('../esther_paper/hybrid_gfdl-esm4_ssp126_default_yield_soybean_min_shift_2015_2098','.tif', sep = "")) +names(raster_2012_hybrid) <- 'hybrid_shift' +# Process the data and convert to csv according to globiom grid +v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) +output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) +globiom_grid_am$shift_hybrid <- output_hybrid +shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) +names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') +shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + +my_parts <- strsplit(test, "/") +paste(my_parts[[1]][1], '/', my_parts[[1]][2], '/output_csv/', my_parts[[1]][3],'.csv', sep = "") diff --git a/code/globiom_harvest_area_converter.R b/code/globiom_harvest_area_converter.R new file mode 100644 index 0000000..bd9da0a --- /dev/null +++ b/code/globiom_harvest_area_converter.R @@ -0,0 +1,42 @@ +library(ncdf4) # package for netcdf manipulation +library(raster) # package for raster manipulation +library(rgdal) # package for geospatial analysis +library(ggplot2) # package for plotting +library(terra) +library(tidyr) +library(dplyr) +library(sf) +library(sp) + +rdir <- ("C:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data/output_shocks_am/scenarios_forum") +setwd(rdir) + +# Load global grid matrix GLOBIOM is based on: +globiom_grid <- st_read("../../COLROW/COLROW_GLOBIOMcode.shp") +# plot(globiom_grid) + +# Load a given globiom csv data (harvest area 2050) +globiom_grid_2050 <- read.csv("c:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/other_papers/Soy_Area_2050.csv") %>% dplyr::select(CR_ID,X2050) %>% mutate(X2050 = X2050 * 1000) %>% rename(area = X2050) + +# merge the global grid with the data above on a common variable, here called 'CR_ID' +globiom_grid_new_area <- merge(globiom_grid, globiom_grid_2050, by='CR_ID') +st_write(globiom_grid_new_area, "soy_harvest_area_globiom_2050.shp") +#plot(globiom_grid_new_area) + +summary(globiom_grid_new_area$area ) + +# Load reference: +refe_grid <- raster('hybrid_ukesm1-0-ll_ssp585_default_yield_soybean_shift_2017-2044.tif') + +ext <- extent(globiom_grid_new_area) +rr <- raster(ext, res=0.5) +crs(rr) <- "+proj=longlat" +rrr <- rasterize(globiom_grid_new_area, rr,"area") + +rrr_resamp <- resample(rrr, refe_grid) + +writeRaster(rrr_resamp,"soy_harvest_area_globiom_2050.tiff",overwrite=TRUE) + +writeRaster(rrr_resamp,"soy_harvest_area_globiom_2050.nc",overwrite=TRUE, format="CDF", varname="area", varunit="ha", + longname="harvest area", xname="lon", yname="lat") + diff --git a/code/globiom_mask_generator.R b/code/globiom_mask_generator.R new file mode 100644 index 0000000..3b18984 --- /dev/null +++ b/code/globiom_mask_generator.R @@ -0,0 +1,50 @@ +library(ncdf4) # package for netcdf manipulation +library(raster) # package for raster manipulation +library(rgdal) # package for geospatial analysis +library(ggplot2) # package for plotting +library(terra) +library(tidyr) +library(dplyr) +library(sf) +library(sp) + +rdir <- ("C:/users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data") +setwd(rdir) + +globiom_grid_am <- vect("./soy_americas_harvest_area.shp") + +# Short code to generate masks for globiom + + +### MAIN CODE +create_shifters_globiom_2012 <- function(path) { + + # Load .nc file + nc_file <- raster(paste(path,'.nc', sep = "")) + + writeRaster(x = nc_file, filename = paste(path,'.tif', sep = ""), format = 'GTiff', overwrite = TRUE) + + raster_2012_hybrid <- rast(paste(path,'.tif', sep = "")) + names(raster_2012_hybrid) <- 'hybrid_shift' + + v_all <- terra::extract(raster_2012_hybrid, globiom_grid_am, fun = mean) + + output_hybrid = unlist(lapply(v_all$hybrid_shift, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) + + globiom_grid_am$shift_hybrid <- output_hybrid + + shift_2012_soybean_am <- data.frame(globiom_grid_am$CR_ID, globiom_grid_am$shift_hybrid) + names(shift_2012_soybean_am) <- c('ALLCOLROW','hybrid') + shift_2012_soybean_am_sub <- na.omit(shift_2012_soybean_am, cols = c("hybrid")) + + + # Save shifters to csv + write.csv(shift_2012_soybean_am_sub, paste(path,'.csv', sep = ""), row.names = FALSE) + +} + +# main code, choose file desired +path_1 = 'soybean_yields_america_detrended_2012' +path_2 = 'soybean_harvested_area_america_2012' + +create_shifters_globiom_2012(path_2) diff --git a/code/gtsm_tests.py b/code/gtsm_tests.py new file mode 100644 index 0000000..862cd95 --- /dev/null +++ b/code/gtsm_tests.py @@ -0,0 +1,107 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Sep 12 17:58:11 2022 + +@author: morenodu +""" + +import os +os.chdir('D:/paper_3/code') +import xarray as xr +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd +import geopandas as gpd +import cartopy +import cartopy.crs as ccrs +import cartopy.io.shapereader as shpreader +import matplotlib as mpl + +mpl.rcParams['axes.spines.right'] = False +mpl.rcParams['axes.spines.top'] = False +mpl.rcParams['figure.dpi'] = 144 +mpl.rcParams.update({'font.size': 14}) + + +# Import data +''' +here we explore the storm surge & tide levels for the Xynthia (2010) and Xaver (2013) stomrs +''' + +DS_storm_surge = xr.open_dataset("../data/surge_tides/reanalysis_waterlevel_hourly_2013_12_v1.nc") +DS_storm_surge_xaver = DS_storm_surge.sel(time = slice('2013-12-05','2013-12-10')) + +DS_storm_surge_xaver['waterlevel'].sel(stations = 2).plot(x='station_x_coordinate', y='station_y_coordinate',transform=ccrs.PlateCarree(), robust=True) + +rm = {'station_x_coordinate':'lon', 'station_y_coordinate':'lat'} + +DS_storm_surge_xaver_test = DS_storm_surge_xaver.rename(rm) +DS_storm_surge_xaver_test['waterlevel'].isel(time = 0, stations = 1).plot(x='lon', y='lat') + +## CODEC +bbox = 44.182204,-6.899414,54.188155,9.887695 +bbox_hamburg = 4.987793,52.133488,13.820801,57.326521 + +rm = {'station_x_coordinate':'lon', 'station_y_coordinate':'lat'} + +ds_gtsm0 = xr.open_dataset('../data/surge_tides/reanalysis_waterlevel_hourly_2013_12_v1.nc').rename(rm) +# Create MultiIndex coordinate +ds_gtsm0_2 = ds_gtsm0.set_index(station=["lat", "lon"]).sel(lat=slice(bbox[1],bbox[3]), lon=slice(bbox[0],bbox[2])) +# Unstack the MultiIndex +ds_gtsm0_2 = ds_gtsm0_2.unstack() + + + +ds_gtsm0.waterlevel.isel(time=0).to_series().plot(x='lon', y='lat') + + + +# Now try with gridded data: +ds_gtsm_gridded = xr.open_dataset('../data/surge_tides/era5-water_level-2013-m12-v0.0.nc').rename(rm) + +# Create MultiIndex coordinate +stations_2 = ds_gtsm_gridded.stations.values +ds_gtsm_gridded = ds_gtsm_gridded.assign_coords(stations = ('stations',stations_2)) +ds_multiindex = ds_gtsm_gridded.set_index(station=["lat", "lon"]).sel(lat=slice(bbox[1],bbox[3]), lon=slice(bbox[0],bbox[2])) +# Unstack the MultiIndex +ds_gtsm_gridded_2 = ds_multiindex.unstack() + +ds_gtsm_gridded_2['waterlevel'].isel(time = 0).plot(x='lon', y='lat',transform=ccrs.PlateCarree(), robust=True) + + +# construct a full grid +def expand_grid(lat,lon): + '''list all combinations of lats and lons using expand_grid(lat,lon)''' + test = [(A,B) for A in lat for B in lon] + test = np.array(test) + test_lat = test[:,0] + test_lon = test[:,1] + full_grid = pd.DataFrame({'lat': test_lat, 'lon': test_lon}) + full_grid = full_grid.sort_values(by=['lat','lon']) + full_grid = full_grid.reset_index(drop=True) + return full_grid + + +data_onto_full_grid = pd.merge(full_grid, your_actual_data,how='left") + + + + + +# compare two datasets: +ds_gtsm0.sel(stations = 3985) +ds_gtsm_gridded.sel(stations = 3985) + + +plt.plot(ds_gtsm0['waterlevel'].sel(stations = 3985), label = '2022_dataset') +plt.plot(ds_gtsm_gridded['waterlevel'].sel(stations = 3990), color = 'red', label = '2020_dataset') +plt.legend() +plt.show() + +plt.figure(figsize = (10,10)) +plt.title('Xaver') +plt.plot(ds_gtsm0['waterlevel'].sel(stations = 3985, time = slice('2013-12-05','2013-12-12')).time,ds_gtsm0['waterlevel'].sel(stations = 3985, time = slice('2013-12-05','2013-12-12')), label = '2022_dataset' ) +plt.plot(ds_gtsm_gridded['waterlevel'].sel(stations = 3990, time = slice('2013-12-05','2013-12-12')).time,ds_gtsm_gridded['waterlevel'].sel(stations = 3990, time = slice('2013-12-05','2013-12-12')) , label = '2020_dataset') +plt.legend() +plt.show() + diff --git a/code/ml_predictions_auto_am.py b/code/ml_predictions_auto_am.py index b419ae5..cd627e2 100644 --- a/code/ml_predictions_auto_am.py +++ b/code/ml_predictions_auto_am.py @@ -493,422 +493,3 @@ def country_location_add(df, US = DS_y_obs_us_all['Yield'], BR = DS_y_obs_br['Yi df_input_hybrid_fut_ipsl_585_am, df_prediction_proj_ipsl_585_am = hybrid_predictions_future(DS_y_epic_proj_ipsl_585_am, df_feature_proj_ipsl_585_am, model = 'ipsl', rcp_scenario = 'ssp585', region = "_am", hybrid_model_full = full_model_hyb_am2, co2_scen = co2_scen) df_input_hybrid_fut_ipsl_126_am, df_prediction_proj_ipsl_126_am = hybrid_predictions_future(DS_y_epic_proj_ipsl_126_am, df_feature_proj_ipsl_126_am, model = 'ipsl', rcp_scenario = 'ssp126', region = "_am", hybrid_model_full = full_model_hyb_am2, co2_scen = co2_scen) - - - - -# DS_y_epic_proj_ukesm_585_am_test_us = DS_y_epic_proj_ukesm_585_am.where(DS_historical_hybrid_us['Yield'].mean('time') > -10) - -# DS_feature_proj_ukesm_585_am = xr.Dataset.from_dataframe(df_feature_proj_ukesm_585_am) -# DS_feature_proj_ukesm_585_am = rearrange_latlot(DS_feature_proj_ukesm_585_am) -# DS_feature_proj_ukesm_585_am_test_us = DS_feature_proj_ukesm_585_am.where(DS_historical_hybrid_us['Yield'].mean('time') > -10) -# df_feature_proj_ukesm_585_am_test_us = DS_feature_proj_ukesm_585_am_test_us.to_dataframe() - - - -# df_input_hybrid_fut_ukesm_585_am, df_prediction_proj_ukesm_585_am = hybrid_predictions_future(DS_y_epic_proj_ukesm_585_am_test_us, df_feature_proj_ukesm_585_am_test_us, model = 'ukesm', rcp_scenario = 'ssp585', region = "_am", hybrid_model_full = full_model_hyb2_am_rf, co2_scen = co2_scen) - - - - - - - - - - - - - - - - - -# #%% FUNCTION TO GENERATE PROJECTIONS BASED ON RANDOM FOREST -# def past_generation_hybrid(region, hybrid_model_full, start_date, end_date, three_models = False, sim_round = '\Gen_Assem'): -# if region == '_am' or 'am': - -# #COMBINE THE THREE REGIONS WITH SOUTH AMERICA 12 MONTHS SHIFTED -# DS_clim_ext_projections_us = xr.open_mfdataset('../../paper_hybrid_agri/data/climpact-master/climpact-master/www/output_historical_us/monthly_data/*.nc').sel(time=slice(start_date, end_date)) - -# # Shift br data one year ahead -# DS_clim_ext_projections_br =xr.open_mfdataset('../../paper_hybrid_agri/data/climpact-master/climpact-master/www/output_gswp3/monthly_data/*.nc').sel(time=slice(start_date, end_date)) - -# DS_clim_ext_projections_br = DS_clim_ext_projections_br.copy().shift(time = 12) # SHIFT EPIC BR ONE YEAR FORWARD - -# # Shift ARG data one year ahead -# DS_clim_ext_projections_arg = xr.open_mfdataset('../../paper_hybrid_agri/data/climpact-master/climpact-master/www/output_historical_arg/monthly_data/*.nc').sel(time=slice(start_date, end_date)) -# DS_clim_ext_projections_arg = DS_clim_ext_projections_arg.copy().shift(time = 12) # SHIFT EPIC BR ONE YEAR FORWARD - -# # Combine all grids -# DS_clim_ext_projections = DS_clim_ext_projections_us.combine_first(DS_clim_ext_projections_br) -# DS_clim_ext_projections = DS_clim_ext_projections.combine_first(DS_clim_ext_projections_arg) -# DS_clim_ext_projections = DS_clim_ext_projections.sel(time=slice(start_date, end_date)) - -# # Reindex to avoid missnig coordinates and dimension values -# DS_clim_ext_projections = rearrange_latlot(DS_clim_ext_projections) - -# #control -# plot_2d_am_map(DS_clim_ext_projections['txm'].isel(time = 0), title = 'First year projection') -# plot_2d_am_map(DS_clim_ext_projections['txm'].isel(time = -1), title = 'Lat year projection') - -# ### Climatic variables -# # Clean -# DS_clim_ext_projections = DS_clim_ext_projections.drop_vars(['fd','id','time_bnds']) # Always zero -# # Selected features -# list_features_am = list_features # From historical time -# DS_clim_ext_projections = DS_clim_ext_projections[list_features_am] - -# DS_clim_ext_projections = DS_clim_ext_projections.where(DS_y_obs_am_det_regul.mean('time') >= -5.0) -# plot_2d_am_map(DS_clim_ext_projections['prcptot'].mean('time')) - -# da_list = [] -# for feature in list(DS_clim_ext_projections.keys()): -# if (type(DS_clim_ext_projections[feature].values[0,0,0]) == np.timedelta64): -# print('Time') -# DS = timedelta_to_int(DS_clim_ext_projections, feature) -# else: -# print('Integer') -# DS = DS_clim_ext_projections[feature] -# da_list.append(DS) - -# DS_clim_ext_projections_combined = xr.merge(da_list) - -# DS_clim_ext_projections_combined = rearrange_latlot(DS_clim_ext_projections_combined) -# DS_clim_ext_projections_combined = DS_clim_ext_projections_combined.reindex(lat=DS_clim_ext_projections_combined.lat[::-1]) -# if len(DS_clim_ext_projections_combined.coords) >3 : -# DS_clim_ext_projections_combined = DS_clim_ext_projections_combined.drop('spatial_ref') - -# DS_clim_ext_projections_am = DS_clim_ext_projections_combined.where( DS_chosen_calendar_am >= 0 ) - -# # ============================================================================= -# # CONVERT CLIMATIC VARIABLES ACCORDING TO THE SOYBEAN GROWING SEASON PER GRIDCELL -# # ============================================================================= -# df_clim_proj_twoyears = dynamic_calendar(DS_clim_ext_projections_am) - -# ### Select specific months -# suffixes = tuple(["_"+str(j) for j in range(1,4)]) -# df_feature_proj_6mon = df_clim_proj_twoyears.loc[:,df_clim_proj_twoyears.columns.str.endswith(suffixes)] - -# df_feature_proj_6mon = df_feature_proj_6mon.rename_axis(index={'year':'time'}) -# df_feature_proj_6mon = df_feature_proj_6mon.reorder_levels(['time','lat','lon']).sort_index() -# df_feature_proj_6mon = df_feature_proj_6mon.dropna() -# df_feature_proj_6mon = df_feature_proj_6mon.astype(float) - -# # ============================================================================= -# # # SECOND DETRENDING PART - SEASONAL CYCLE -# # # If to keep historical mean values: DS_feature_season_6mon_am_det, if to keep projections mean: DS_feature_proj_6mon_am -# # ============================================================================= -# DS_feature_proj_6mon_am = xr.Dataset.from_dataframe(df_feature_proj_6mon) -# DS_feature_proj_6mon_am = rearrange_latlot(DS_feature_proj_6mon_am) - -# # Choose historical mean so there is no mean deviation between the historical and future timelines -# DS_feature_proj_6mon_am_det = detrend_dataset(DS_feature_proj_6mon_am, deg = 'free', mean_data = DS_feature_season_6mon_am_det.sel(time = slice(2000,2015))) - -# df_feature_proj_6mon_am_det = DS_feature_proj_6mon_am_det.to_dataframe().dropna() -# df_feature_proj_6mon_am_det = df_feature_proj_6mon_am_det.rename_axis(list(DS_feature_proj_6mon_am_det.coords)).reorder_levels(['time','lat','lon']).sort_index() - -# for feature in df_feature_proj_6mon.columns: -# df_feature_proj_6mon[feature].groupby('time').mean().plot(label = 'old') -# df_feature_proj_6mon_am_det[feature].groupby('time').mean().plot(label = 'detrend') -# plt.title(f'{feature} for past') -# plt.legend() -# plt.show() - -# list_feat_precipitation = [s for s in df_feature_proj_6mon_am_det.keys() if "prcptot" in s] -# for feature in list_feat_precipitation: -# df_feature_proj_6mon_am_det[feature][df_feature_proj_6mon_am_det[feature] < 0] = 0 - -# # UPDATE DETRENDED VALUES -# df_feature_proj_6mon = df_feature_proj_6mon_am_det -# # DS_features = xr.Dataset.from_dataframe(df_feature_proj_6mon) -# # DS_features.to_netcdf('output_clim_indices_'+region +'/clim_indices_'+model_full+'_'+rcp_scenario+'_soybean_2015_2100.nc') - -# # ============================================================================= -# # #%% EPIC projections -# # ============================================================================= -# def epic_projections_function_co2(co2_scenario): -# if region == '_am' or 'am': #Write down a way of selecting each file and merging them toegether with the right windows - -# DS_y_epic_proj = xr.open_dataset("epic-iiasa_gswp3-w5e5_obsclim_2015soc_default_yield-soy-noirr_global_annual_1901_2016.n", decode_times=False) -# DS_y_epic_proj = convert_timeunit(DS_y_epic_proj) - -# # EPIC DATA USA - keep as it is -# DS_y_epic_proj_us = DS_y_epic_proj.where(DS_y_epic_us['yield'].mean('time') >= -5.0 ) - -# # Shift br data one year ahead -# DS_y_epic_proj_br = DS_y_epic_proj.where(DS_y_epic_br['yield'].mean('time') >= -5.0 ) -# DS_y_epic_proj_br = DS_y_epic_proj_br.shift(time = 1) # SHIFT EPIC BR ONE YEAR FORWARD - -# # Shift ARG data one year ahead -# DS_y_epic_proj_arg = DS_y_epic_proj.where(DS_y_epic_arg['yield'].mean('time') >= -5.0 ) -# DS_y_epic_proj_arg = DS_y_epic_proj_arg.shift(time = 1) # SHIFT EPIC BR ONE YEAR FORWARD - -# # Combine all grids -# DS_y_epic_proj_am = DS_y_epic_proj_us.combine_first(DS_y_epic_proj_br) -# DS_y_epic_proj_am = DS_y_epic_proj_am.combine_first(DS_y_epic_proj_arg) - -# # Reindex to avoid missnig coordinates and dimension values -# DS_y_epic_proj_am = DS_y_epic_proj_am.reindex(lat=DS_y_epic_proj_am.lat[::-1]) -# DS_y_epic_proj_am = DS_y_epic_proj_am.sortby(['time','lat','lon']) -# DS_y_epic_proj_am = rearrange_latlot(DS_y_epic_proj_am) - -# DS_y_epic_proj_am = DS_y_epic_proj_am.where(DS_y_epic_am_det_regul.mean('time') >= -5.0 ) - -# else: -# DS_y_epic_proj = xr.open_dataset("epic-iiasa_past_w5e5_yield-soy-noirr_global_annual_1900_1970.nc", decode_times=False) -# DS_y_epic_proj_am = DS_y_epic_proj.where(DS_y_epic_am_det_regul.mean('time') >= -5.0 ) - -# return DS_y_epic_proj_am, DS_feature_proj_6mon_am, df_feature_proj_6mon -# # ============================================================================= - -# DS_y_epic_proj_am, DS_feature_proj_6mon_am, df_feature_proj_6mon2 = epic_projections_function_co2(co2_scenario = co2_scen) - - -# return DS_y_epic_proj_am, DS_feature_proj_6mon_am, df_feature_proj_6mon2 - -# # Prediction past-pre training - -# DS_y_epic_proj_ukesm_585_am, DS_feature_proj_ukesm_585_am, df_feature_proj_ukesm_585_am = past_generation_hybrid(region = "_am", hybrid_model_full = model_to_be_used, start_date='1900-12-12', end_date='1972-12-12') - - -# # #%% -# # # UKESM model -# # df_input_hybrid_fut_ukesm_585_am, df_prediction_proj_ukesm_585_am = projections_generation_hybrid(model = 'ukesm', rcp_scenario = 'ssp585', region = "_am", hybrid_model_full = model_to_be_used, start_date='01-01-2015', end_date='31-12-2100', co2_scen = co2_scen) -# # df_input_hybrid_fut_ukesm_126_am, df_prediction_proj_ukesm_126_am = projections_generation_hybrid(model = 'ukesm', rcp_scenario = 'ssp126', region = "_am", hybrid_model_full = model_to_be_used, start_date='01-01-2015', end_date='31-12-2100', co2_scen = co2_scen) - -# # # GFDL model -# # df_input_hybrid_fut_gfdl_585_am, df_prediction_proj_gfdl_585_am = projections_generation_hybrid(model = 'gfdl', rcp_scenario = 'ssp585', region = "_am", hybrid_model_full = model_to_be_used, start_date='01-01-2015', end_date='31-12-2100', co2_scen = co2_scen) -# # df_input_hybrid_fut_gfdl_126_am, df_prediction_proj_gfdl_126_am = projections_generation_hybrid(model = 'gfdl', rcp_scenario = 'ssp126', region = "_am", hybrid_model_full = model_to_be_used, start_date='01-01-2015', end_date='31-12-2100', co2_scen = co2_scen) - -# # # IPSL model -# # df_input_hybrid_fut_ipsl_585_am, df_prediction_proj_ipsl_585_am = projections_generation_hybrid(model = 'ipsl', rcp_scenario = 'ssp585', region = "_am", hybrid_model_full = model_to_be_used, start_date='01-01-2015', end_date='31-12-2100', co2_scen = co2_scen) -# # df_input_hybrid_fut_ipsl_126_am, df_prediction_proj_ipsl_126_am = projections_generation_hybrid(model = 'ipsl', rcp_scenario = 'ssp126', region = "_am", hybrid_model_full = model_to_be_used, start_date='01-01-2015', end_date='31-12-2100', co2_scen = co2_scen) - - -# # #%% TEST partial dependence plot for out of calibration zone cases - 2088 should be extreme - -# # load future prodictions -# df_predict_fut = df_prediction_proj_ukesm_585_am.iloc[:,[0]].copy() # Use 2015 case -# df_predict_fut = df_predict_fut.rename(columns={'yield-soy-noirr_default':'yield-soy-noirr'}) - -# df_proj_fut = df_input_hybrid_fut_ukesm_585_am.rename(columns={'yield-soy-noirr_default':'yield-soy-noirr'}) -# df_proj_fut = df_proj_fut[ ['yield-soy-noirr'] + [ col for col in df_proj_fut.columns if col != 'yield-soy-noirr' ] ] -# # df_proj_fut = df_proj_fut.drop(columns='yield-soy-noirr_default') - -# df_hybrid_am_test = df_input_hybrid_am.copy() -# df_hybrid_am_test = df_hybrid_am_test.rename(columns={'yield':'yield-soy-noirr'}) - -# # EPIC -# plt.plot(df_hybrid_am_test['yield-soy-noirr'].groupby('time').mean(), label = 'History') -# plt.plot(df_proj_fut['yield-soy-noirr'].groupby('time').mean(), label = 'Future') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().idxmin(), linestyle = 'dashed') -# plt.title('EPIC predictions') -# plt.legend() -# plt.show() - -# # TMX -# plt.plot(df_hybrid_am_test['txm_4'].groupby('time').mean(), label = 'History') -# plt.plot(df_proj_fut['txm_4'].groupby('time').mean(), label = 'Future') -# plt.title('TXM_4 predictions') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().idxmin(), linestyle = 'dashed') -# plt.axvline(df_proj_fut['txm_4'].groupby('time').mean().idxmax(), c = 'r', linestyle = 'dashed') -# plt.legend() -# plt.show() - -# # HYBRID -# plt.plot(df_predict_hyb_am['Yield'].groupby('time').mean(), label = 'History') -# plt.plot(df_predict_fut['yield-soy-noirr'].groupby('time').mean(), label = 'Future') -# plt.title('Hybrid predictions') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().idxmin(), linestyle = 'dashed') -# plt.legend() -# plt.show() - -# # Plots for points distribution -# for feature in df_proj_fut.columns: -# df_clim_extrapolated = df_proj_fut[feature].where(df_proj_fut[feature] > df_hybrid_am_test[feature].max()).dropna() -# df_y_extrapolated = df_predict_fut['yield-soy-noirr'].where(df_proj_fut[feature] > df_hybrid_am_test[feature].max()).dropna() - -# plt.scatter(df_hybrid_am_test[feature], df_predict_hyb_am['Yield'], color = 'k', label = 'History') -# plt.scatter(df_proj_fut[feature], df_predict_fut['yield-soy-noirr'], alpha = 0.8, label = 'Projection') -# sns.regplot(df_hybrid_am_test[feature], df_predict_hyb_am['Yield'], color = 'k', label = 'History', scatter = False) -# sns.regplot(df_proj_fut[feature], df_predict_fut['yield-soy-noirr'], label = 'Projection', scatter = False) -# plt.scatter(df_clim_extrapolated, df_y_extrapolated, alpha = 0.8, label = 'Extrapolation') -# plt.legend(loc="upper right") -# plt.title(f'Scatterplot of {feature} for GCM-RCPs') -# plt.ylabel('Yield') -# if feature in ['txm_3','txm_4','txm_5']: -# x_label = 'Temperature (°C)' -# elif feature in ['prcptot_2', 'prcptot_3', 'prcptot_4', 'prcptot_5']: -# x_label = 'Precipitation (mm/month)' -# else: -# x_label = 'Yield (ton/ha)' - -# plt.xlabel(x_label) -# plt.show() - -# for feature in df_proj_fut.columns: -# sns.kdeplot(df_hybrid_am_test[feature],fill=True, alpha = 0.3, label = 'History') -# sns.kdeplot(df_proj_fut[feature],fill=True, alpha = 0.3, label = 'Proj') -# print('hist', np.round(df_hybrid_am_test[feature].mean(), 3), 'fut', np.round(df_proj_fut[feature].mean(),3)) -# plt.legend() -# plt.show() - -# plt.plot(df_hybrid_am_test[feature].groupby('time').mean(), label = 'History') -# plt.plot(df_proj_fut[feature].groupby('time').mean(), label = 'Future') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().idxmin(), linestyle = 'dashed') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().nsmallest(5).index[1], linestyle = 'dashed') -# plt.axvline(df_predict_fut['yield-soy-noirr'].groupby('time').mean().idxmin(), color = 'red',linestyle = 'dashed') -# plt.title(f'{feature} predictions') -# plt.legend() -# plt.show() - -# sns.kdeplot(df_predict_hyb_am['Yield'], fill=True, alpha = 0.3, label = 'History') -# sns.kdeplot(df_predict_fut['yield-soy-noirr'],fill=True, alpha = 0.3, label = 'Proj') -# plt.title('Hybrid predictions') -# plt.legend() -# plt.show() - -# plt.plot(df_predict_hyb_am['Yield'].groupby('time').mean(), label = 'History') -# plt.plot(df_predict_fut['yield-soy-noirr'].groupby('time').mean(), label = 'Future') -# plt.title('Hybrid predictions') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().idxmin(), linestyle = 'dashed') -# plt.axvline(df_proj_fut['yield-soy-noirr'].groupby('time').mean().nsmallest(5).index[1], linestyle = 'dashed') -# plt.legend() -# plt.show() - -# # for feature in df_proj_fut.columns: -# # df_clim_extrapolated = df_proj_fut[feature].where(df_proj_fut[feature] < df_hybrid_am_test[feature].min()).dropna() -# # df_y_extrapolated = df_predict_fut['yield-soy-noirr'].where(df_proj_fut[feature] < df_hybrid_am_test[feature].min()).dropna() - -# # plt.scatter(df_hybrid_am_test[feature], df_predict_hyb_am['Yield'], color = 'k') -# # plt.scatter(df_proj_fut[feature], df_predict_fut['yield-soy-noirr'], alpha = 0.8) -# # plt.hlines(df_predict_hyb_am['Yield'].mean(), df_hybrid_am_test[feature].min(), df_hybrid_am_test[feature].max(), color = 'k') -# # plt.scatter(df_clim_extrapolated, df_y_extrapolated, alpha = 0.8) -# # # plt.legend(loc="upper right") -# # plt.title(f'Scatterplot of {feature} for GCM-RCPs') -# # plt.ylabel('Yield') -# # if feature in ['tnx_3','tnx_4','tnx_5']: -# # x_label = 'Temperature (°C)' -# # elif feature in ['prcptot_3', 'prcptot_4', 'prcptot_5']: -# # x_label = 'Precipitation (mm/month)' -# # else: -# # x_label = 'Yield (ton/ha)' - -# # plt.xlabel(x_label) -# # plt.show() - - -# #%% Partial dependence plots -# from sklearn.inspection import PartialDependenceDisplay - -# features_to_plot = ['yield-soy-noirr','txm_3'] -# fig, ax1 = plt.subplots(1, 1, figsize=(10, 8), dpi=500) -# disp1 = PartialDependenceDisplay.from_estimator(model_to_be_used, df_proj_fut, features_to_plot, pd_line_kw={'color':'r'},percentiles=(0.01,0.99), ax = ax1) -# disp2 = PartialDependenceDisplay.from_estimator(model_to_be_used, df_hybrid_am_test, features_to_plot, ax = disp1.axes_,percentiles=(0,1), pd_line_kw={'color':'k'}) -# plt.setp(disp1.deciles_vlines_, visible=False) -# plt.setp(disp2.deciles_vlines_, visible=False) -# plt.show() - -# fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 8)) -# disp1.plot(ax=[ax1, ax2], line_kw={"label": "Extrapolation", "color": "red"}) -# disp2.plot(ax=[ax1, ax2], line_kw={"label": "Training", "color": "black"}) -# ax1.set_ylim(1.0, 2.8) -# ax2.set_ylim(1.0, 2.8) -# plt.setp(disp1.deciles_vlines_, visible=False) -# plt.setp(disp2.deciles_vlines_, visible=False) -# ax1.legend() -# plt.show() - -# features_to_plot = [1,2,3] -# fig, ax1 = plt.subplots(1, 1, figsize=(10, 8), dpi=500) -# disp3 = PartialDependenceDisplay.from_estimator(model_to_be_used, df_proj_fut, features_to_plot, pd_line_kw={'color':'r'},percentiles=(0.01,0.99), ax = ax1) -# disp4 = PartialDependenceDisplay.from_estimator(model_to_be_used, df_hybrid_am_test, features_to_plot, ax = disp3.axes_,percentiles=(0,1), pd_line_kw={'color':'k'}) -# plt.ylim(0, 2.6) -# plt.setp(disp3.deciles_vlines_, visible=False) -# plt.setp(disp4.deciles_vlines_, visible=False) - -# fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8)) -# disp3.plot(ax=[ax1, ax2, ax3], line_kw={"label": "Extrapolation", "color": "red"}) -# disp4.plot( -# ax=[ax1, ax2, ax3], line_kw={"label": "Training", "color": "black"} -# ) -# ax1.set_ylim(1, 2.8) -# ax2.set_ylim(1, 2.8) -# ax3.set_ylim(1, 2.8) -# plt.setp(disp3.deciles_vlines_, visible=False) -# plt.setp(disp4.deciles_vlines_, visible=False) -# ax1.legend() -# plt.show() - -# features_to_plot = [4,5,6] -# fig, ax1 = plt.subplots(1, 1, figsize=(10, 8), dpi=500) -# disp5 = PartialDependenceDisplay.from_estimator(model_to_be_used, df_proj_fut, features_to_plot, pd_line_kw={'color':'r'},percentiles=(0.01,0.99), ax = ax1, method = 'brute') -# disp6 = PartialDependenceDisplay.from_estimator(model_to_be_used, df_hybrid_am_test, features_to_plot, ax = disp5.axes_,percentiles=(0,1), pd_line_kw={'color':'k'}, method = 'brute') -# plt.ylim(0, 2.6) -# plt.setp(disp5.deciles_vlines_, visible=False) -# plt.setp(disp6.deciles_vlines_, visible=False) - -# fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8)) -# disp5.plot(ax=[ax1, ax2, ax3], line_kw={"label": "Extrapolation", "color": "red"}) -# disp6.plot( -# ax=[ax1, ax2, ax3], line_kw={"label": "Training", "color": "black"} -# ) -# ax1.set_ylim(1, 2.8) -# ax2.set_ylim(1, 2.8) -# ax3.set_ylim(1, 2.8) -# plt.setp(disp5.deciles_vlines_, visible=False) -# plt.setp(disp6.deciles_vlines_, visible=False) -# ax1.legend() -# plt.show() - -# #%% -# mean_stuff = df_proj_fut['prcptot_4'].groupby(['lat','lon']).mean().to_xarray() -# stuff_2088 = df_proj_fut['prcptot_4'].loc[2088].groupby(['lat','lon']).mean().to_xarray() -# delta = stuff_2088 - mean_stuff -# delta.plot() -# # test_values = df_hybrid_am_test.iloc[[0],:-1].copy() -# # test_values.iloc[[0],:] = [0,0,0,0,50,50,50] -# # test_prdict = full_model_hyb.predict(test_values) -# # print(test_prdict) - -# DS_test = xr.open_mfdataset('monthly_ukesm_ssp585_am/txm_MON_climpact.ukesm1-0-ll_r1i1p1f2_w5e5_ssp5-ssp5.nc').where(DS_y_epic_am_det_regul.mean('time') >= -5.0) - -# # txm 2 is peaking one year after txm_4 -# DS_test['txm'].mean(['lat','lon']).plot() -# DS_test['txm'].sel(time = slice('2030-01-16','2040-10-16')).mean(['lat','lon']).plot() -# DS_test['txm'].sel(time = slice('2037-01-01','2037-01-30')).mean(['lat','lon']).values -# DS_test['txm'].sel(time = slice('2036-11-01','2037-03-30')).mean(['lat','lon']).values - -# DS_test['txm'].sel(time = slice('2034-11-01','2035-03-30')).mean(['lat','lon']).values -# DS_test['txm'].sel(time = slice('2035-11-01','2036-03-30')).mean(['lat','lon']).values # March SHould be equal to TXM_4 peak -# DS_test['txm'].sel(time = slice('2036-11-01','2037-03-30')).mean(['lat','lon']).values # January SHould be equal to TXM_2 peak - -# df_test = DS_test['txm'].mean(['lat','lon']).to_dataframe() -# df_test - -# from statsmodels.tsa.seasonal import seasonal_decompose -# decompose_data = seasonal_decompose(df_test, model="additive", period = 516) -# decompose_data.plot(); - -# decompose_data.seasonal.plot(); - -# from statsmodels.graphics.tsaplots import plot_acf -# plot_acf(df_test); - -# plt.figure(figsize=(10,6)) -# plt.show() - -# DS_test['txm'].sel(time = '2100-12-16').mean(['time']).plot() - - -# DS_test2 = xr.open_mfdataset('epic-iiasa_ukesm1-0-ll_w5e5_ssp585_2015soc_2015co2_yield-soy-noirr_global_annual_2015_2100.nc', decode_times= False) -# DS_test2 = DS_test2.where(DS_y_epic_am_det.mean('time') >= -5.0) - - -# DS_test2['yield-soy-noirr'].sel(time = 438).plot() -# plot_2d_am_map(DS_test2['yield-soy-noirr'].sel(time = 438)) -# #%% - - - diff --git a/code/paper2_results_comparison_production.py b/code/paper2_results_comparison_production.py index 6d54f5b..249531b 100644 --- a/code/paper2_results_comparison_production.py +++ b/code/paper2_results_comparison_production.py @@ -491,6 +491,7 @@ def convert_clim_weighted_ensemble(df_clim, DS_counterfactuals_weighted_country, ax.hlines(y=df_climatology_grouped[name], xmin=df_climatology_grouped['month']-1.5, xmax=df_climatology_grouped['month']-0.5, colors='firebrick', linestyles='--', lw=2, label = '2012 event') g1 = sns.boxplot(ax = ax, x="month", y=name, data=df_analogues_grouped, palette = ['tab:orange'], hue = 'scenario') #drawstyle = 'steps-mid', err_kws= {'step':'mid'} + g1 = sns.swarmplot(ax = ax, x="month", y=name, data=df_analogues_grouped, palette = ['firebrick'], size = 10, alpha=0.5, hue = 'scenario') #drawstyle = 'steps-mid', err_kws= {'step':'mid'} # g1.set(xticklabels=[]) # remove the tick labels # g1.set(xlabel= '') g1.set(ylabel= '' ) @@ -504,7 +505,7 @@ def convert_clim_weighted_ensemble(df_clim, DS_counterfactuals_weighted_country, axes[1].set_title("b) Temperature",loc='left') handles, labels = ax.get_legend_handles_labels() -fig.legend(handles, labels, bbox_to_anchor=(0.8, 0.01), ncol=3, frameon=False) +fig.legend(handles[:2], labels[:2], bbox_to_anchor=(0.8, 0.01), ncol=3, frameon=False) plt.tight_layout() plt.savefig('paper_figures_production/clim_conditions_grouped_2012.png', format='png', dpi=300, bbox_inches='tight') plt.show() @@ -1393,6 +1394,7 @@ def clim_conditions_analogues_anom_2(DS_area_hist, DS_area_fut, DS_counterfactua ax.hlines(y=data_country_climatology.query('year == 2012')[name], xmin=data_country_climatology.query('year == 2012')['month']-1.5, xmax=data_country_climatology.query('year == 2012')['month']-0.5, colors='firebrick', linestyles='--', lw=2, label = '2012 event') g1 = sns.boxplot(ax = ax, x="month", y=name, data=data_country, hue = 'level', palette=['tab:orange']) #drawstyle = 'steps-mid', err_kws= {'step':'mid'} + # g1 = sns.swarmplot(ax = ax, x="month", y=name, data=data_country, hue = 'level', palette=['firebrick']) #drawstyle = 'steps-mid', err_kws= {'step':'mid'} g1.set(ylabel= '' ) ax.get_legend().remove() diff --git a/code/paper_2_americas_sim2.py b/code/paper_2_americas_sim2.py index bb5debc..d11fb9c 100644 --- a/code/paper_2_americas_sim2.py +++ b/code/paper_2_americas_sim2.py @@ -19,6 +19,8 @@ import cartopy.crs as ccrs import cartopy.io.shapereader as shpreader import matplotlib as mpl +import rasterio +from rasterio.features import geometry_mask mpl.rcParams['axes.spines.right'] = False mpl.rcParams['axes.spines.top'] = False @@ -186,28 +188,63 @@ def weighted_conversion(DS, DS_area, name_ds = 'Yield'): DS_weighted = ((DS * DS_area['harvest_area'].where(DS > -10) / DS_area['harvest_area'].where(DS > -10).sum(['lat','lon']))) return DS_weighted.sum(['lat','lon']) + +# List of countries +countries = shpreader.natural_earth(resolution='50m',category='cultural',name='admin_0_countries') +gdf = gpd.read_file(countries) +import odc.geo.xr + +def country_sel(da, country = 'USA'): + if type(country) == str: + country = [country] + gdf_country = gdf[gdf['ADM0_A3'].isin(country)] + + # Convert the polygon to a shapely object + polygon = gdf_country.geometry.unary_union + + # Create a mask of the polygon bounds + mask = geometry_mask([polygon], transform=da.odc.geobox.transform, + out_shape=da.odc.geobox.shape, + all_touched=False, invert=True) + + # Use the mask to extract the data within the bounds of the polygon + da_subset = da.where(mask) + return da_subset + #%% LOADING MAIN DATA # Load Country shapes -us1_shapes = states_mask('../../Paper_drought/data/gadm36_USA_1.shp') -br1_shapes = states_mask('../../Paper_drought/data/gadm36_BRA_1.shp') -arg_shapes = states_mask('GIS/gadm36_ARG_1.shp') +us1_shapes = states_mask('../../Paper_hybrid_agri/data/countries_shapefiles/gadm36_USA_1.shp') +br1_shapes = states_mask('../../Paper_hybrid_agri/data/countries_shapefiles/gadm36_BRA_1.shp') +arg_shapes = states_mask('../../Paper_hybrid_agri/data/countries_shapefiles/gadm36_ARG_1.shp') # ============================================================================= -# USe MIRCA to isolate the rainfed 90% soybeans + # Use SPAM2010 to isolate the soybean grid cells where rainfed soybean is > 0 # ============================================================================= # DS_mirca_test = xr.open_dataset("../../paper_hybrid_agri/data/americas_mask_ha.nc", decode_times=False).rename({'latitude': 'lat', 'longitude': 'lon','annual_area_harvested_rfc_crop08_ha_30mn':'harvest_area'}) DS_mirca_test = xr.open_dataset("../../paper_hybrid_agri/data/soy_harvest_spam_native_05x05.nc", decode_times=False) plot_2d_am_map(DS_mirca_test['harvest_area']) #### HARVEST DATA -DS_harvest_area_sim = xr.load_dataset("../../paper_hybrid_agri/data/soybean_harvest_area_calculated_americas_hg.nc", decode_times=False) -DS_harvest_area_sim = DS_harvest_area_sim.sel(time=2012) -DS_harvest_area_sim = DS_harvest_area_sim.where(DS_mirca_test['harvest_area'] > 0 ) +DS_harvest_area_total = xr.load_dataset("../../paper_hybrid_agri/data/soybean_harvest_area_calculated_americas_hg.nc", decode_times=False) +DS_harvest_area_total2012 = DS_harvest_area_total.sel(time=2012) +plot_2d_am_map(DS_harvest_area_total2012['harvest_area']) + +DS_harvest_area_sim = DS_harvest_area_total2012.where(DS_mirca_test['harvest_area'] > 0 ) DS_harvest_area_sim = rearrange_latlot(DS_harvest_area_sim) DS_mirca_test = DS_harvest_area_sim plot_2d_am_map(DS_mirca_test['harvest_area']) # plot_2d_am_map(DS_harvest_area_sim['harvest_area']) +print('fraction of rainfed areas from total:', DS_harvest_area_sim['harvest_area'].sum().values/DS_harvest_area_total2012['harvest_area'].sum().values) + +def rainfed_ratio(DS_y_mask): + test_rf = DS_harvest_area_sim['harvest_area'].where(DS_y_mask['Yield'] > -10) + test_total = DS_harvest_area_total2012['harvest_area'].where(DS_y_mask['Yield'] > -10) + print('fraction of rainfed areas from total:', test_rf.sum().values/test_total.sum().values) + + plot_2d_am_map(test_rf) + plot_2d_am_map(test_total) + # ============================================================================= # GLOBAL - EPIC # ============================================================================= @@ -219,7 +256,8 @@ def weighted_conversion(DS, DS_area, name_ds = 'Yield'): DS_y_epic = DS_y_epic.sel(time=slice('1972-12-12','2016-12-12')) DS_y_epic = DS_y_epic.rename({'yield-soy-noirr':'yield'}) -DS_y_epic_am = DS_y_epic.where(DS_mirca_test['harvest_area'] > 0 ) +DS_y_epic_am = country_sel(DS_y_epic, country = ['USA','BRA','ARG']) +# DS_y_epic_am = DS_y_epic.where(DS_mirca_test['harvest_area'] > 0 ) plot_2d_am_map(DS_y_epic_am['yield'].isel(time=0)) plot_2d_am_map(DS_y_epic_am['yield'].sel(time=2016)) # remove weird data point on north of brazil that seems to be off-calendar @@ -228,6 +266,12 @@ def weighted_conversion(DS, DS_area, name_ds = 'Yield'): DS_y_epic_am['yield'] = DS_y_epic_am['yield'].transpose('time','lat','lon') DS_y_epic_am = rearrange_latlot(DS_y_epic_am) +DS_y_epic_us = country_sel(DS_y_epic_am, country = 'USA') +DS_y_epic_br = country_sel(DS_y_epic_am, country = 'BRA') +DS_y_epic_arg = country_sel(DS_y_epic_am, country = 'ARG') + + + # ============================================================================= # OBSERVED CENSUS DATA # ============================================================================= @@ -237,44 +281,53 @@ def weighted_conversion(DS, DS_area, name_ds = 'Yield'): # Convert time unit --- # units, reference_date = DS_y_obs_us_all.time.attrs['units'].split('since') #pd.date_range(start=reference_date, periods=DS_y_obs_us_all.sizes['time'], freq='YS').year DS_y_obs_us_all['time'] = DS_y_obs_us_all['time'].astype(int) DS_y_obs_us_all = DS_y_obs_us_all.sel(time=slice(start_date, end_date)) +DS_y_obs_us_2012 = DS_y_obs_us_all.sel(time = 2012) DS_y_obs_us_all = DS_y_obs_us_all.where(DS_mirca_test['harvest_area'] > 0 ) DS_y_obs_us_all.to_netcdf("soybean_yields_US_1978_2016.nc") plot_2d_am_map(DS_y_obs_us_all['Yield'].mean('time')) -DS_y_epic_us = DS_y_epic_am.where(DS_y_obs_us_all['Yield'] > - 5) +# DS_y_epic_us = DS_y_epic_am.where(DS_y_obs_us_all['Yield'] > - 5) # ============================================================================= # # BRAZIL -------------------------------------------------------- # ============================================================================= DS_y_obs_br = xr.open_dataset("../../paper_hybrid_agri/data/soy_yield_1975_2016_05x05_1prc.nc", decode_times=False) DS_y_obs_br=DS_y_obs_br.sel(time = slice(start_date, end_date)) +DS_y_obs_br_2012 = DS_y_obs_br.sel(time = 2012) DS_y_obs_br = DS_y_obs_br.where(DS_mirca_test['harvest_area'] > 0 ) plot_2d_am_map(DS_y_obs_br['Yield'].mean('time')) DS_y_obs_br.to_netcdf("soybean_yields_BR_1978_2016.nc") # SHIFT EPIC FOR BRAZIL ONE YEAR FORWARD TO MATCH INTERNATIONAL CALENDARS -DS_y_epic_br = DS_y_epic_am.where(DS_y_obs_br['Yield'].mean('time') > - 5) +# DS_y_epic_br = DS_y_epic_am.where(DS_y_obs_br['Yield'].mean('time') > - 5) DS_y_epic_br = DS_y_epic_br.copy().shift(time = 1) # SHIFT EPIC BR ONE YEAR FORWARD -DS_y_epic_br = DS_y_epic_br.where(DS_y_obs_br['Yield'] > - 5) +# DS_y_epic_br = DS_y_epic_br.where(DS_y_obs_br['Yield'] > - 5) # plot_2d_am_map(DS_y_epic_br['yield'].sel(time=1979)) # ============================================================================= # # AREGNTINA -------------------------------------------------------- # ============================================================================= DS_y_obs_arg = xr.open_dataset("../../paper_hybrid_agri/data/soy_yield_arg_1974_2019_05x05.nc", decode_times=False)#soy_yield_1980_2016_1prc05x05 / soy_yield_1980_2016_all_filters05x05 -DS_y_obs_arg=DS_y_obs_arg.sel(time = slice(start_date-1, end_date)) # get one year earlier to shift them forward +DS_y_obs_arg = DS_y_obs_arg.sel(time = slice(start_date-1, end_date)) # get one year earlier to shift them forward # SHIFT OBSERVED DATA FOR ARGENTINA ONE YEAR FORWARD TO MATCH INTERNATIONAL CALENDARS - from planting year to harvest year DS_y_obs_arg = DS_y_obs_arg.copy().shift(time = 1) # SHIFT AGRNEITNA ONE YeAR FORWARD +DS_y_obs_arg_2012 = DS_y_obs_arg.sel(time = 2012) + DS_y_obs_arg = DS_y_obs_arg.where(DS_mirca_test['harvest_area'] > 0 ).sel(time=slice(start_date, end_date)) DS_y_obs_arg.to_netcdf("soybean_yields_ARG_1978_2016.nc") # SHIFT EPIC DATA FOR ARGENTINA ONE YEAR FORWARD TO MATCH INTERNATIONAL CALENDARS -DS_y_epic_arg = DS_y_epic_am.where(DS_y_obs_arg['Yield'].mean('time') > - 5) +# DS_y_epic_arg = DS_y_epic_am.where(DS_y_obs_arg['Yield'].mean('time') > - 5) DS_y_epic_arg = DS_y_epic_arg.copy().shift(time = 1) # SHIFT EPIC ARG ONE YeAR FORWARD -DS_y_epic_arg = DS_y_epic_arg.where(DS_y_obs_arg['Yield'] > - 5) +# DS_y_epic_arg = DS_y_epic_arg.where(DS_y_obs_arg['Yield'] > - 5) + +# Calculate the amount of rainfed regions left from the total +rainfed_ratio(DS_y_obs_us_2012) +rainfed_ratio(DS_y_obs_br_2012) +rainfed_ratio(DS_y_obs_arg_2012) # ============================================================================= # # Plots for analysis @@ -333,7 +386,12 @@ def weighted_conversion(DS, DS_area, name_ds = 'Yield'): # Detrend timeseries DS_y_obs_am_det = xr.DataArray( detrend_dim(DS_y_obs_am_clip['Yield'], 'time') + DS_y_obs_am_clip['Yield'].mean('time'), name= DS_y_obs_am_clip['Yield'].name, attrs = DS_y_obs_am_clip['Yield'].attrs) DS_y_obs_am_det = DS_y_obs_am_det.sel(time = slice(start_date_det, end_date_det)) +# Save files DS_y_obs_am_det.to_netcdf("soybean_yields_america_detrended_1978_2016.nc") +DS_y_obs_am_det.sel(time = 2012).to_netcdf("soybean_yields_america_detrended_2012.nc") + +DS_harvest_area_2012 = DS_harvest_area_sim.where(DS_y_obs_am_det.sel(time = 2012) > -10) +DS_harvest_area_2012.to_netcdf("soybean_harvested_area_america_2012.nc") DS_y_epic_am_det = xr.DataArray( detrend_dim(DS_y_epic_am_clip["yield"], 'time') + DS_y_epic_am_clip["yield"].mean('time'), name= DS_y_epic_am_clip["yield"].name, attrs = DS_y_epic_am_clip["yield"].attrs) DS_y_epic_am_det = DS_y_epic_am_det.sel(time = slice(start_date_det, end_date_det)) @@ -439,12 +497,7 @@ def calibration(X_origin,y_origin,type_of_model='RF', test_year = None): ]) elif type_of_model == 'DNN': -# extra layer batch epoch nodes lr dropout_value best_epoch R2 MAE RMSE -# 31 True 256 700 512 0.01 0.2 514 0.651 0.24641 0.34128 -# 1 False 1024 700 512 0.01 0.2 564 0.648 0.24794 0.34267 -# 30 True 256 700 512 0.005 0.2 635 0.648 0.24502 0.34266 -# number of epochs 529 or 431 - epochs_train = 441 #390 + epochs_train = 441 batch_size_train = 1024 nodes_size = 512 learning_rate_train = 0.01 @@ -486,10 +539,6 @@ def create_model(): train_model.compile(loss='mse', optimizer=tf.keras.optimizers.Adam(learning_rate_train), metrics=['mean_squared_error','mean_absolute_error']) return train_model - # Callbacks to monitor the performance of the optimization of the model and if there is any overfitting - # callback_model = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience= 100, restore_best_weights=True) - # mc = ModelCheckpoint('best_model_test.h5', monitor='val_loss', mode='min', save_best_only=True, verbose=1) - model_rf = Pipeline([ ('scaler', StandardScaler()), ('estimator', KerasRegressor(model=create_model(), epochs= epochs_train, random_state = 0, batch_size=batch_size_train, verbose=0)) ]) #, callbacks=callback_model #, validation_split= 0.1, callbacks=[callback_model,mc] @@ -606,6 +655,12 @@ def units_conversion(DS_exclim): # COMBINE DS_exclim_am = DS_exclim_us.combine_first(DS_exclim_arg) DS_exclim_am = DS_exclim_am.combine_first(DS_exclim_br) + +# FOR ESTHER +ds_harvest_area_2050 = xr.open_dataset("output_shocks_am/scenarios_forum/soy_harvest_area_globiom_2050.nc", decode_times=False) +ds_harvest_area_2050 = ds_harvest_area_2050.isel(lat=slice(None, None, -1)) +DS_exclim_am_mask = DS_exclim_am.where(ds_harvest_area_2050['area'] > -10) + DS_exclim_am = DS_exclim_am.where(DS_y_obs_am_det.mean('time') > -10) DS_exclim_am = rearrange_latlot(DS_exclim_am) @@ -819,7 +874,7 @@ def dynamic_calendar(DS, df_calendar_month_am): print('Dynamic ECE results:') df_feature_season_6mon_am_countryloc = country_location_add(df_feature_season_6mon_am) X, y = df_feature_season_6mon_am_countryloc, df_obs_am_det_clip['Yield'] -y_pred_exclim_dyn_am, y_pred_total_exclim_dyn_am, model_exclim_dyn_am, full_model_exclim_dyn_am = calibration(X, y, type_of_model='DNN') +y_pred_exclim_dyn_am, y_pred_total_exclim_dyn_am, model_exclim_dyn_am, full_model_exclim_dyn_am = calibration(X, y, type_of_model='RF') y_pred_exclim_dyn_am_2012, y_pred_total_exclim_dyn_am_2012, model_exclim_dyn_am_2012, full_model_exclim_dyn_am_2012 = calibration(X, y, type_of_model='RF',test_year = 2012) # # ============================================================================= @@ -875,7 +930,7 @@ def dynamic_calendar(DS, df_calendar_month_am): # Standard model print('Standard Epic results:') -y_pred_epic_am, y_pred_total_epic_am, model_epic_am, full_model_epic_am = calibration(X, y, type_of_model='DNN') +y_pred_epic_am, y_pred_total_epic_am, model_epic_am, full_model_epic_am = calibration(X, y, type_of_model='RF') y_pred_epic_am_am_2012, y_pred_total_epic_am_2012, model_epic_am_2012, full_model_epic_am_2012 = calibration(X, y, type_of_model='RF',test_year = 2012) #%% Hybrid model @@ -921,7 +976,7 @@ def dynamic_calendar(DS, df_calendar_month_am): # Deep neural network - takes some time to run y_pred_hyb_am2, y_pred_total_hyb_am2, model_hyb_am2, full_model_hyb_am2 = calibration(X, y, type_of_model='DNN') -y_pred_hyb_am_2012, y_pred_total_hyb_am_2012, model_hyb_am_2012, full_model_hyb_am_2012 = calibration(X, y, type_of_model='RF',test_year = 2012) +# y_pred_hyb_am_2012, y_pred_total_hyb_am_2012, model_hyb_am_2012, full_model_hyb_am_2012 = calibration(X, y, type_of_model='RF',test_year = 2012) # # ============================================================================= # # Save the Model diff --git a/code/soy_harvest_spam_globiom.csv b/code/soy_harvest_spam_globiom.csv new file mode 100644 index 0000000..883026d --- /dev/null +++ b/code/soy_harvest_spam_globiom.csv @@ -0,0 +1,2381 @@ +"ALLCOLROW","harvest_area" +"CR235206",1.60000002384186 +"CR237206",261.186553955078 +"CR237207",77.8345794677734 +"CR238205",237.481018066406 +"CR238206",594.807067871094 +"CR238207",381.575775146484 +"CR239205",161.233581542969 +"CR239206",273.640716552734 +"CR239207",772.716369628906 +"CR239208",261.423767089844 +"CR240207",460.069427490234 +"CR240208",323.872955322266 +"CR240209",196.839599609375 +"CR240210",231.523422241211 +"CR240211",432.213409423828 +"CR240212",150.973937988281 +"CR240213",378.59814453125 +"CR241208",568.979248046875 +"CR241209",246.103652954102 +"CR241210",224.448104858398 +"CR241211",299.788970947266 +"CR241212",61.3138389587402 +"CR241213",190.789733886719 +"CR242210",51.5463638305664 +"CR242211",196.421661376953 +"CR242212",78.8677062988281 +"CR242213",189.716354370117 +"CR243206",2976.2861328125 +"CR243207",3725.53247070312 +"CR243210",175.227798461914 +"CR243211",1821.74328613281 +"CR243212",4781.05078125 +"CR243213",104.572601318359 +"CR244202",4.00075435638428 +"CR244203",4.09846448898315 +"CR244204",1038.43640136719 +"CR244205",929.246154785156 +"CR244206",1552.75183105469 +"CR244207",2694.2216796875 +"CR244210",234.450927734375 +"CR244211",606.21044921875 +"CR244212",152.671112060547 +"CR244213",48.846809387207 +"CR244214",10.0894193649292 +"CR244220",273.320037841797 +"CR245201",3.79999995231628 +"CR245202",5.20470237731934 +"CR245203",16.7932147979736 +"CR245204",17.759105682373 +"CR245205",626.807739257812 +"CR245209",2253.71801757812 +"CR245210",671.992370605469 +"CR245211",45.5729331970215 +"CR245212",39.8870124816895 +"CR245213",17.0595226287842 +"CR245214",8.34681987762451 +"CR245217",23.2999992370605 +"CR245218",382.827697753906 +"CR245219",281.250946044922 +"CR245241",74.4549865722656 +"CR246202",2.40000009536743 +"CR246203",276.175842285156 +"CR246204",4.19999980926514 +"CR246205",809.744812011719 +"CR246206",937.123718261719 +"CR246207",1051.39721679688 +"CR246209",1396.51037597656 +"CR246210",417.740814208984 +"CR246211",97.9846725463867 +"CR246212",30.1667175292969 +"CR246213",24.7593441009521 +"CR246214",1.79999995231628 +"CR246240",2.86467623710632 +"CR246241",1.39999997615814 +"CR247200",15.6999998092651 +"CR247201",12.1439113616943 +"CR247202",172.199996948242 +"CR247203",1086.8525390625 +"CR247204",246.09553527832 +"CR247205",1641.12585449219 +"CR247206",1832.86120605469 +"CR247207",1349.208984375 +"CR247208",1740.13171386719 +"CR247209",1660.14221191406 +"CR247210",247.217163085938 +"CR247211",315.645629882812 +"CR247212",62.4353370666504 +"CR247213",110.78116607666 +"CR247214",4.26068878173828 +"CR247215",14.8999996185303 +"CR247221",5.74941968917847 +"CR247222",146.705551147461 +"CR247223",265.39404296875 +"CR247224",199.936294555664 +"CR247239",245.878662109375 +"CR247240",2.35844802856445 +"CR247241",11.5106630325317 +"CR248199",4.1112961769104 +"CR248200",12.3198146820068 +"CR248201",162.382995605469 +"CR248202",155.910186767578 +"CR248203",1038.35266113281 +"CR248204",1075.2314453125 +"CR248205",2477.67407226562 +"CR248206",4280.7587890625 +"CR248207",4116.11376953125 +"CR248208",2637.68579101562 +"CR248209",2610.15795898438 +"CR248210",664.299438476562 +"CR248211",158.15983581543 +"CR248212",57.5129470825195 +"CR248213",40.5371551513672 +"CR248214",3.42409634590149 +"CR248215",25.8775672912598 +"CR248221",11.3336124420166 +"CR248222",125.032402038574 +"CR248223",164.130462646484 +"CR248224",190.58708190918 +"CR248225",481.822692871094 +"CR248238",270.063751220703 +"CR248239",381.765594482422 +"CR248240",92.7595748901367 +"CR248241",110.8193359375 +"CR248242",173.741790771484 +"CR249200",12.1782293319702 +"CR249201",199.174301147461 +"CR249202",223.547988891602 +"CR249203",916.741271972656 +"CR249204",2639.94213867188 +"CR249205",5002.83447265625 +"CR249206",5304.4296875 +"CR249207",4479.98583984375 +"CR249208",3335.66333007812 +"CR249209",1036.54321289062 +"CR249210",475.283660888672 +"CR249211",237.477554321289 +"CR249212",136.244781494141 +"CR249213",30.4386329650879 +"CR249214",26.754524230957 +"CR249215",42.556957244873 +"CR249221",42.2214050292969 +"CR249222",33.6214904785156 +"CR249225",1803.18603515625 +"CR249237",553.836669921875 +"CR249238",674.321716308594 +"CR249239",688.304809570312 +"CR249240",200.745330810547 +"CR249241",254.367904663086 +"CR249242",208.487899780273 +"CR249243",312.899993896484 +"CR250200",275.143890380859 +"CR250201",696.380493164062 +"CR250202",496.817504882812 +"CR250203",517.882202148438 +"CR250204",1554.88500976562 +"CR250205",2659.26953125 +"CR250206",2655.70874023438 +"CR250207",3762.46044921875 +"CR250208",1807.89929199219 +"CR250209",837.888549804688 +"CR250210",84.7434997558594 +"CR250211",1207.83715820312 +"CR250212",1861.71557617188 +"CR250213",777.319274902344 +"CR250214",94.1346817016602 +"CR250215",202.081130981445 +"CR250216",417.783355712891 +"CR250217",46.9467658996582 +"CR250219",66.8610305786133 +"CR250220",63.3186073303223 +"CR250221",78.4744567871094 +"CR250236",166.417297363281 +"CR250237",1353.25329589844 +"CR250238",686.363647460938 +"CR250239",797.152893066406 +"CR250240",402.621307373047 +"CR250241",536.374267578125 +"CR250242",211.379898071289 +"CR250243",253.367706298828 +"CR251200",66.1704788208008 +"CR251201",538.304321289062 +"CR251202",460.484619140625 +"CR251203",567.404846191406 +"CR251204",787.979614257812 +"CR251205",745.587158203125 +"CR251206",829.324462890625 +"CR251207",1695.85754394531 +"CR251208",1521.97888183594 +"CR251209",378.654418945312 +"CR251210",249.035354614258 +"CR251211",2543.931640625 +"CR251212",854.470458984375 +"CR251213",1296.53308105469 +"CR251214",1221.06848144531 +"CR251215",2291.7529296875 +"CR251216",1417.59240722656 +"CR251217",198.150360107422 +"CR251218",129.680969238281 +"CR251219",1180.34436035156 +"CR251220",754.397399902344 +"CR251221",91.5736999511719 +"CR251222",867.146789550781 +"CR251232",4270.466796875 +"CR251236",2052.9970703125 +"CR251237",3186.34619140625 +"CR251238",2360.103515625 +"CR251239",1257.40173339844 +"CR251240",525.458068847656 +"CR251241",699.649658203125 +"CR251242",540.990600585938 +"CR251243",333.757019042969 +"CR252200",46.373893737793 +"CR252201",58.4483680725098 +"CR252202",28.8184547424316 +"CR252203",192.077865600586 +"CR252204",375.887878417969 +"CR252205",746.274353027344 +"CR252206",142.589797973633 +"CR252207",299.014312744141 +"CR252208",282.916687011719 +"CR252209",504.288116455078 +"CR252212",960.313537597656 +"CR252213",662.80908203125 +"CR252214",1125.6865234375 +"CR252215",2244.4765625 +"CR252216",809.302978515625 +"CR252217",222.699813842773 +"CR252218",784.027404785156 +"CR252220",580.262634277344 +"CR252221",153.934310913086 +"CR252222",113.082489013672 +"CR252229",3980.15625 +"CR252230",2437.34716796875 +"CR252231",3227.93920898438 +"CR252232",394.007751464844 +"CR252235",760.611633300781 +"CR252236",3730.88549804688 +"CR252237",4253.78076171875 +"CR252238",4259.5361328125 +"CR252239",2864.25659179688 +"CR252240",558.036499023438 +"CR252241",1176.48596191406 +"CR252242",338.430267333984 +"CR252243",198.310424804688 +"CR252244",240.114562988281 +"CR253200",37.392993927002 +"CR253201",59.3216857910156 +"CR253202",37.2000007629395 +"CR253203",433.916198730469 +"CR253207",445.103302001953 +"CR253208",429.147064208984 +"CR253209",453.757049560547 +"CR253212",938.698486328125 +"CR253213",558.502868652344 +"CR253214",1818.43603515625 +"CR253215",1082.06018066406 +"CR253216",198.108276367188 +"CR253220",66.7712554931641 +"CR253222",25.8373470306396 +"CR253226",385.516143798828 +"CR253227",434.241119384766 +"CR253228",506.731506347656 +"CR253229",4269.3681640625 +"CR253230",4301.84228515625 +"CR253231",2147.82495117188 +"CR253232",2157.8388671875 +"CR253233",565.521240234375 +"CR253234",204.829467773438 +"CR253235",1394.04211425781 +"CR253236",4331.3603515625 +"CR253237",4605.21875 +"CR253238",4624.77734375 +"CR253239",3484.02783203125 +"CR253240",1139.18896484375 +"CR253241",839.329833984375 +"CR253242",210.107116699219 +"CR253243",114.290107727051 +"CR253244",155.807907104492 +"CR253245",214.902755737305 +"CR254207",950.429077148438 +"CR254208",620.952392578125 +"CR254209",232.001602172852 +"CR254210",1322.17065429688 +"CR254211",1490.8046875 +"CR254212",606.382995605469 +"CR254213",804.494262695312 +"CR254214",389.573028564453 +"CR254215",383.672149658203 +"CR254216",1362.69177246094 +"CR254219",404.345642089844 +"CR254220",120.441772460938 +"CR254222",16.0747871398926 +"CR254223",19.4379768371582 +"CR254226",264.696136474609 +"CR254227",708.96142578125 +"CR254228",571.588317871094 +"CR254229",4177.09375 +"CR254230",4936.01025390625 +"CR254231",2436.90234375 +"CR254232",2104.82983398438 +"CR254233",1638.72717285156 +"CR254234",624.100646972656 +"CR254235",916.004211425781 +"CR254236",3561.46118164062 +"CR254237",4244.22021484375 +"CR254238",3995.49658203125 +"CR254239",1975.44750976562 +"CR254240",874.830200195312 +"CR254241",1032.47326660156 +"CR254242",402.874114990234 +"CR254243",213.312149047852 +"CR254244",95.8080825805664 +"CR254245",337.691497802734 +"CR254246",206.048252741496 +"CR254247",18.4554271697998 +"CR254248",35.7361450195312 +"CR255201",241.657440185547 +"CR255202",288.415649414062 +"CR255203",200.810760498047 +"CR255204",744.965087890625 +"CR255205",1189.14086914062 +"CR255206",1074.06848144531 +"CR255207",913.683959960938 +"CR255208",509.305908203125 +"CR255209",427.809631347656 +"CR255210",948.945007324219 +"CR255211",393.972747802734 +"CR255212",333.281768798828 +"CR255213",818.034240722656 +"CR255214",517.32958984375 +"CR255215",1411.44641113281 +"CR255216",956.514099121094 +"CR255219",549.385864257812 +"CR255220",114.95499420166 +"CR255222",33.3204193115234 +"CR255223",1.38181102275848 +"CR255226",1417.85583496094 +"CR255227",462.741119384766 +"CR255228",1143.16357421875 +"CR255229",5211.92724609375 +"CR255230",2215.4716796875 +"CR255231",1152.28088378906 +"CR255232",2770.41479492188 +"CR255233",3038.39453125 +"CR255234",1034.17321777344 +"CR255235",1590.56079101562 +"CR255236",3657.97534179688 +"CR255237",4066.83422851562 +"CR255238",3714.7734375 +"CR255239",722.553833007812 +"CR255240",890.317138671875 +"CR255241",738.662292480469 +"CR255242",107.632530212402 +"CR255243",305.044708251953 +"CR255244",303.019989013672 +"CR255245",389.460052490234 +"CR255246",75.9972410202026 +"CR255247",13.506552696228 +"CR256186",2.20356965065002 +"CR256188",1.20000004768372 +"CR256200",551.812927246094 +"CR256201",252.649978637695 +"CR256202",360.469635009766 +"CR256203",517.193420410156 +"CR256204",876.429016113281 +"CR256205",1893.2646484375 +"CR256206",1819.72778320312 +"CR256207",1263.14416503906 +"CR256208",785.938110351562 +"CR256209",446.072784423828 +"CR256210",534.144348144531 +"CR256211",387.67724609375 +"CR256212",342.794097900391 +"CR256213",70.6897735595703 +"CR256214",363.1201171875 +"CR256215",1214.52880859375 +"CR256216",1902.00329589844 +"CR256220",70.6669235229492 +"CR256221",70.7426071166992 +"CR256222",71.44580078125 +"CR256223",3.69044351577759 +"CR256224",63.8319664001465 +"CR256226",283.9580078125 +"CR256227",1961.80151367188 +"CR256228",3913.40576171875 +"CR256229",3273.26098632812 +"CR256230",1703.07092285156 +"CR256231",1415.8427734375 +"CR256232",2248.02490234375 +"CR256233",2648.18725585938 +"CR256234",1990.93615722656 +"CR256235",119.407371520996 +"CR256236",2790.1826171875 +"CR256237",4184.00537109375 +"CR256238",1545.87463378906 +"CR256239",235.891006469727 +"CR256240",532.653991699219 +"CR256241",518.374755859375 +"CR256242",77.41796875 +"CR256243",329.406616210938 +"CR256244",316.948455810547 +"CR256245",43.8194847106934 +"CR257200",682.053955078125 +"CR257201",376.547485351562 +"CR257203",2311.7587890625 +"CR257204",656.428161621094 +"CR257205",1059.81604003906 +"CR257206",302.94189453125 +"CR257207",299.967102050781 +"CR257208",298.315368652344 +"CR257209",91.9593658447266 +"CR257210",467.231872558594 +"CR257211",24.1000003814697 +"CR257212",240.601593017578 +"CR257213",126.490158081055 +"CR257214",392.791931152344 +"CR257215",1278.18103027344 +"CR257219",99.0737533569336 +"CR257220",11.307412147522 +"CR257221",11.0889987945557 +"CR257222",57.6249656677246 +"CR257223",98.1409443855286 +"CR257224",86.1528625488281 +"CR257226",567.952026367188 +"CR257227",2654.79833984375 +"CR257228",2510.65991210938 +"CR257229",2138.03857421875 +"CR257230",2236.62939453125 +"CR257231",2050.91455078125 +"CR257232",1830.85046386719 +"CR257233",550.359069824219 +"CR257234",664.118896484375 +"CR257235",357.199340820312 +"CR257236",2416.71533203125 +"CR257237",2424.49584960938 +"CR257238",562.563354492188 +"CR257239",553.749084472656 +"CR257240",506.300354003906 +"CR257241",315.185852050781 +"CR257242",181.54670715332 +"CR257243",224.776107788086 +"CR258200",596.028015136719 +"CR258201",1040.20007324219 +"CR258202",152.664443969727 +"CR258204",587.511596679688 +"CR258205",512.772766113281 +"CR258206",60.5416069030762 +"CR258207",5.91680288314819 +"CR258208",9.36883163452148 +"CR258209",28.0672588348389 +"CR258210",29.3684902191162 +"CR258211",28.5209369659424 +"CR258212",222.046752929688 +"CR258213",114.411109924316 +"CR258214",371.681030273438 +"CR258215",3393.78979492188 +"CR258219",14.9591684341431 +"CR258220",7.55647039413452 +"CR258221",21.6432647705078 +"CR258222",143.339477539062 +"CR258223",446.153137207031 +"CR258226",2139.81323242188 +"CR258227",3552.8095703125 +"CR258228",2510.34619140625 +"CR258229",1179.68005371094 +"CR258230",965.749145507812 +"CR258231",802.2275390625 +"CR258232",265.030883789062 +"CR258233",213.554336547852 +"CR258234",124.270027160645 +"CR258238",882.032836914062 +"CR258239",1502.55004882812 +"CR258240",11.1508617401123 +"CR259200",171.572006225586 +"CR259201",531.921325683594 +"CR259202",110.149993896484 +"CR259203",87.182243347168 +"CR259204",238.698348999023 +"CR259205",20.1000003814697 +"CR259207",3.74305701255798 +"CR259208",3.87987661361694 +"CR259209",7.88086032867432 +"CR259211",56.5508346557617 +"CR259212",100.168731689453 +"CR259214",896.628967285156 +"CR259217",1197.31115722656 +"CR259218",162.24040222168 +"CR259221",88.6061325073242 +"CR259222",105.735359191895 +"CR259223",117.159721374512 +"CR259224",78.5333938598633 +"CR259225",410.297973632812 +"CR259226",2460.47534179688 +"CR259227",3526.88305664062 +"CR259228",2175.04028320312 +"CR259229",1103.56774902344 +"CR259230",1556.544921875 +"CR259231",1539.51379394531 +"CR259232",1369.71032714844 +"CR259233",1516.31506347656 +"CR259239",751.99462890625 +"CR259240",511.119903564453 +"CR259241",31.16259765625 +"CR260201",34.4000015258789 +"CR260206",10.1141147613525 +"CR260207",18.9703903198242 +"CR260208",24.0305919647217 +"CR260212",51.8539123535156 +"CR260213",101.586242675781 +"CR260214",473.648864746094 +"CR260215",1100.27673339844 +"CR260216",1953.97155761719 +"CR260220",137.000350952148 +"CR260221",48.8162612915039 +"CR260222",124.110237121582 +"CR260223",390.430938720703 +"CR260224",26.984130859375 +"CR260225",299.470062255859 +"CR260226",2128.56079101562 +"CR260227",2002.20373535156 +"CR260228",921.246459960938 +"CR260229",2313.71118164062 +"CR260230",2773.65625 +"CR260231",2662.23950195312 +"CR260232",1255.38586425781 +"CR260233",1093.50708007812 +"CR260238",398.746124267578 +"CR261194",9.44968795776367 +"CR261200",273.517456054688 +"CR261201",145.044204711914 +"CR261202",486.468627929688 +"CR261203",94.0624771118164 +"CR261204",334.174102783203 +"CR261205",341.656433105469 +"CR261206",9.86913585662842 +"CR261207",43.6170692443848 +"CR261208",75.7447052001953 +"CR261209",35.996826171875 +"CR261210",15.0805168151855 +"CR261212",236.626724243164 +"CR261213",370.754180908203 +"CR261214",366.058898925781 +"CR261215",1195.75537109375 +"CR261220",50.2730522155762 +"CR261221",109.99242401123 +"CR261222",108.703559875488 +"CR261223",115.602149963379 +"CR261224",31.1122665405273 +"CR261226",665.083557128906 +"CR261227",355.357696533203 +"CR261228",1093.87561035156 +"CR261229",1274.75927734375 +"CR261230",2266.98779296875 +"CR261231",1043.90405273438 +"CR261232",1396.88708496094 +"CR261233",768.497741699219 +"CR261234",190.478805541992 +"CR261235",230.956802368164 +"CR262194",8.29749965667725 +"CR262195",28.6651916503906 +"CR262196",36.5206985473633 +"CR262197",38.6428413391113 +"CR262198",275.934661865234 +"CR262199",386.600006103516 +"CR262202",16.9170875549316 +"CR262203",22.1128311157227 +"CR262204",206.387268066406 +"CR262205",325.389007568359 +"CR262206",85.5900192260742 +"CR262207",71.8675689697266 +"CR262208",99.5420913696289 +"CR262209",929.119079589844 +"CR262210",533.850769042969 +"CR262211",244.19221496582 +"CR262212",47.579833984375 +"CR262213",166.204223632812 +"CR262214",209.483459472656 +"CR262215",781.244934082031 +"CR262216",865.490356445312 +"CR262218",1274.90454101562 +"CR262220",137.313110351562 +"CR262221",165.704162597656 +"CR262222",57.5067710876465 +"CR262223",40.9689140319824 +"CR262224",92.969841003418 +"CR262226",313.101806640625 +"CR262229",796.66064453125 +"CR262230",40.4503135681152 +"CR262231",667.979125976562 +"CR262232",768.4853515625 +"CR262233",289.790496826172 +"CR263196",41.3627738952637 +"CR263197",166.444976806641 +"CR263198",534.599060058594 +"CR263199",256.2841796875 +"CR263200",217.8359375 +"CR263201",380.421722412109 +"CR263202",285.711029052734 +"CR263203",228.547775268555 +"CR263204",76.1044998168945 +"CR263205",240.992553710938 +"CR263206",66.3442001342773 +"CR263207",97.1832885742188 +"CR263208",59.6329345703125 +"CR263209",304.720123291016 +"CR263210",266.871490478516 +"CR263211",451.559234619141 +"CR263212",196.258071899414 +"CR263213",933.546203613281 +"CR263214",1571.06848144531 +"CR263215",876.038696289062 +"CR263216",809.705505371094 +"CR263217",520.233642578125 +"CR263220",663.443359375 +"CR263221",524.795837402344 +"CR263222",212.048187255859 +"CR263224",28.1139698028564 +"CR263225",69.7744064331055 +"CR263226",211.414489746094 +"CR264190",111.991752624512 +"CR264191",1012.16473388672 +"CR264195",37.5999984741211 +"CR264196",435.6064453125 +"CR264197",331.206451416016 +"CR264198",463.692474365234 +"CR264199",879.783020019531 +"CR264200",145.738082885742 +"CR264201",172.156112670898 +"CR264202",470.328796386719 +"CR264203",763.0361328125 +"CR264204",170.114730834961 +"CR264205",77.0117874145508 +"CR264206",66.7668075561523 +"CR264207",68.3734436035156 +"CR264209",148.896713256836 +"CR264210",222.781295776367 +"CR264211",622.68505859375 +"CR264212",425.487762451172 +"CR264213",717.776123046875 +"CR264214",1616.04992675781 +"CR264215",669.8017578125 +"CR264216",1075.92749023438 +"CR264217",259.860778808594 +"CR264220",1127.21240234375 +"CR264221",1592.38916015625 +"CR264222",478.162994384766 +"CR264223",442.861907958984 +"CR264224",83.5796661376953 +"CR264225",23.7031669616699 +"CR265189",224.429565429688 +"CR265190",43.8747444152832 +"CR265191",496.164733886719 +"CR265193",89.3854141235352 +"CR265194",267.46044921875 +"CR265195",142.608139038086 +"CR265196",255.785385131836 +"CR265197",258.335601806641 +"CR265198",176.839324951172 +"CR265199",367.380706787109 +"CR265200",251.431777954102 +"CR265201",742.043151855469 +"CR265202",56.4961891174316 +"CR265203",477.787994384766 +"CR265204",207.010955810547 +"CR265205",127.972671508789 +"CR265206",48.527530670166 +"CR265207",47.153148651123 +"CR265209",267.444702148438 +"CR265210",797.4150390625 +"CR265211",482.888061523438 +"CR265212",629.772705078125 +"CR265213",1401.75512695312 +"CR265214",1804.73718261719 +"CR265215",1792.34228515625 +"CR265216",1299.22729492188 +"CR265217",863.247863769531 +"CR265218",685.052062988281 +"CR265220",1171.97509765625 +"CR265221",662.386352539062 +"CR265222",414.662811279297 +"CR265223",203.348556518555 +"CR265224",38.9057083129883 +"CR265225",39.1384391784668 +"CR266188",223.469024658203 +"CR266189",173.18864440918 +"CR266190",2.93048334121704 +"CR266191",3.18395113945007 +"CR266194",27.8033790588379 +"CR266195",226.865219116211 +"CR266196",196.306503295898 +"CR266198",143.817657470703 +"CR266199",143.44287109375 +"CR266200",12.1420841217041 +"CR266201",6.09999990463257 +"CR266203",26.6569995880127 +"CR266204",291.788208007812 +"CR266205",8.29970550537109 +"CR266206",125.959915161133 +"CR266207",21.1000003814697 +"CR266209",405.151977539062 +"CR266210",313.704162597656 +"CR266211",160.014938354492 +"CR266212",1108.77954101562 +"CR266213",1566.40930175781 +"CR266214",1766.76220703125 +"CR266215",1058.06677246094 +"CR266216",1484.44128417969 +"CR266217",763.743103027344 +"CR266219",658.776000976562 +"CR266220",702.901000976562 +"CR266221",262.253814697266 +"CR266222",93.1477203369141 +"CR266223",50.6226043701172 +"CR266224",161.323196411133 +"CR266225",87.2897262573242 +"CR267187",258.473052978516 +"CR267190",3.10014581680298 +"CR267194",229.806991577148 +"CR267195",433.281188964844 +"CR267196",556.410949707031 +"CR267197",1393.74963378906 +"CR267198",793.250549316406 +"CR267199",521.906127929688 +"CR267200",384.113037109375 +"CR267201",1033.82763671875 +"CR267202",770.426513671875 +"CR267203",638.722961425781 +"CR267204",696.968322753906 +"CR267205",2.68173909187317 +"CR267207",7.71086168289185 +"CR267208",9.65196228027344 +"CR267209",118.158020019531 +"CR267210",232.635589599609 +"CR267211",1081.5078125 +"CR267212",1253.84594726562 +"CR267213",798.288024902344 +"CR267214",937.138671875 +"CR267215",419.975830078125 +"CR267216",622.090454101562 +"CR267217",249.540252685547 +"CR267218",362.399993896484 +"CR267219",374.841247558594 +"CR267220",298.198547363281 +"CR267221",123.824577331543 +"CR267222",54.2893028259277 +"CR267223",81.5366134643555 +"CR267224",44.9000015258789 +"CR267225",99.1727142333984 +"CR268192",93.9738922119141 +"CR268193",76.6769561767578 +"CR268194",429.067169189453 +"CR268195",786.474426269531 +"CR268196",1045.78234863281 +"CR268197",906.624877929688 +"CR268198",976.539306640625 +"CR268199",199.292098999023 +"CR268200",259.3115234375 +"CR268201",524.079467773438 +"CR268202",1362.89172363281 +"CR268203",2018.51904296875 +"CR268204",1920.99975585938 +"CR268205",2602.96533203125 +"CR268206",2765.35034179688 +"CR268207",2431.06127929688 +"CR268208",1061.61975097656 +"CR268209",345.698486328125 +"CR268210",349.776184082031 +"CR268211",660.906188964844 +"CR268212",886.193298339844 +"CR268213",740.59375 +"CR268214",705.313720703125 +"CR268215",482.600891113281 +"CR268216",232.98486328125 +"CR268217",414.235015869141 +"CR268218",188.690963745117 +"CR268219",508.762329101562 +"CR268220",177.562927246094 +"CR268221",53.1915588378906 +"CR268222",57.9683799743652 +"CR268223",60.2704467773438 +"CR268224",11.3000001907349 +"CR269193",141.434234619141 +"CR269194",782.604797363281 +"CR269195",1129.7197265625 +"CR269196",1302.24035644531 +"CR269197",1502.64428710938 +"CR269198",636.4921875 +"CR269199",603.901733398438 +"CR269200",350.343139648438 +"CR269201",1176.20141601562 +"CR269202",1120.95324707031 +"CR269203",1575.91247558594 +"CR269204",1939.30078125 +"CR269205",1957.88061523438 +"CR269206",1846.16259765625 +"CR269207",1275.662109375 +"CR269208",607.550048828125 +"CR269209",399.379577636719 +"CR269210",277.344604492188 +"CR269211",1162.63525390625 +"CR269212",231.360687255859 +"CR269213",120.87287902832 +"CR269215",32.807430267334 +"CR269216",464.192687988281 +"CR269217",312.245178222656 +"CR269218",175.637512207031 +"CR269219",178.746063232422 +"CR269220",56.6668510437012 +"CR269221",198.780075073242 +"CR269222",50.591007232666 +"CR269223",55.5203590393066 +"CR269224",119.287338256836 +"CR270193",143.984222412109 +"CR270194",731.195556640625 +"CR270195",549.887817382812 +"CR270196",1172.03625488281 +"CR270197",1586.67749023438 +"CR270198",992.479370117188 +"CR270199",598.143005371094 +"CR270200",397.003601074219 +"CR270201",93.4653244018555 +"CR270202",245.342514038086 +"CR270203",469.791839599609 +"CR270204",587.996215820312 +"CR270205",948.55078125 +"CR270206",510.965972900391 +"CR270207",420.914093017578 +"CR270208",420.126434326172 +"CR270209",365.981719970703 +"CR270210",230.987915039062 +"CR270211",820.700012207031 +"CR270212",1908.23718261719 +"CR270213",128.63639831543 +"CR270214",368.750183105469 +"CR270215",159.692108154297 +"CR270216",188.770706176758 +"CR270217",118.789039611816 +"CR270218",12.5740480422974 +"CR270220",148.146179199219 +"CR270221",58.968318939209 +"CR270222",73.4000015258789 +"CR270223",78.2611846923828 +"CR270224",131.499130249023 +"CR271193",95.8976287841797 +"CR271194",614.833312988281 +"CR271195",364.537902832031 +"CR271196",1140.72534179688 +"CR271197",1569.97265625 +"CR271198",908.428955078125 +"CR271199",537.727966308594 +"CR271200",314.925476074219 +"CR271201",16.6389484405518 +"CR271202",39.6401901245117 +"CR271204",659.851379394531 +"CR271205",455.013488769531 +"CR271206",512.789794921875 +"CR271207",1729.78454589844 +"CR271208",1126.02734375 +"CR271209",436.402374267578 +"CR271214",421.178680419922 +"CR271215",332.060882568359 +"CR271217",17.2000007629395 +"CR271218",2.29999995231628 +"CR271223",21.4172763824463 +"CR271224",80.6868133544922 +"CR272194",234.511459350586 +"CR272195",674.136657714844 +"CR272196",904.32666015625 +"CR272197",719.4892578125 +"CR272198",939.120422363281 +"CR272199",600.932189941406 +"CR272200",33.3479843139648 +"CR272205",325.166412353516 +"CR272206",659.741760253906 +"CR272207",2067.1787109375 +"CR272208",2018.78845214844 +"CR272209",1407.03698730469 +"CR272217",72.1842651367188 +"CR272218",112.300003051758 +"CR272222",65.3823165893555 +"CR272223",112.750068664551 +"CR272224",95.1655502319336 +"CR273190",22.3810806274414 +"CR273194",144.566253662109 +"CR273195",417.725311279297 +"CR273196",396.700012207031 +"CR273197",185.798721313477 +"CR273206",2138.2431640625 +"CR273207",434.102813720703 +"CR273208",74.4374313354492 +"CR273217",93.7373886108398 +"CR273218",420.561370849609 +"CR273222",78.6999969482422 +"CR274187",186.675582885742 +"CR274188",480.419281005859 +"CR274190",13.2000665664673 +"CR274191",8.31300640106201 +"CR274207",39.0937614440918 +"CR275187",274.568237304688 +"CR275188",850.22119140625 +"CR275189",263.862426757812 +"CR275190",8.66579532623291 +"CR275193",98.859977722168 +"CR275221",128.5 +"CR276187",327.619079589844 +"CR276193",83.168342590332 +"CR228245",45.3866691589355 +"CR228246",187.973907470703 +"CR228247",728.151794433594 +"CR229233",3.52909207344055 +"CR229234",28.5144214630127 +"CR229235",99.9443740844727 +"CR229236",1901.16748046875 +"CR230230",239.798858642578 +"CR230232",564.42529296875 +"CR230233",231.990280151367 +"CR230234",1156.83642578125 +"CR230235",1359.48596191406 +"CR230236",1590.78259277344 +"CR230237",782.912292480469 +"CR230238",703.761108398438 +"CR230245",245.400863647461 +"CR230246",2112.2587890625 +"CR230247",1437.97387695312 +"CR230248",1268.90454101562 +"CR230249",1272.67431640625 +"CR230251",58.1322708129883 +"CR230252",28.0289573669434 +"CR230253",66.2995300292969 +"CR230254",1.69971442222595 +"CR231229",139.480819702148 +"CR231230",146.764190673828 +"CR231231",747.31689453125 +"CR231232",1338.38696289062 +"CR231233",1304.11340332031 +"CR231234",3372.84545898438 +"CR231235",1760.82348632812 +"CR231236",187.579299926758 +"CR231238",124.118530273438 +"CR231239",36.1339721679688 +"CR231241",492.382537841797 +"CR231242",420.878204345703 +"CR231244",1294.55297851562 +"CR231245",733.485290527344 +"CR231246",2438.435546875 +"CR231247",2676.26684570312 +"CR231248",2541.98583984375 +"CR231249",2399.15747070312 +"CR231251",535.818969726562 +"CR231252",492.025390625 +"CR231253",29.4072589874268 +"CR231254",2.3648509979248 +"CR231255",7.28545761108398 +"CR232228",221.694381713867 +"CR232229",481.496459960938 +"CR232230",2281.11010742188 +"CR232231",1036.25805664062 +"CR232232",844.075866699219 +"CR232233",1336.85949707031 +"CR232234",1200.14318847656 +"CR232235",317.359283447266 +"CR232236",123.23070526123 +"CR232237",159.269348144531 +"CR232240",331.488128662109 +"CR232241",483.071228027344 +"CR232242",2006.86096191406 +"CR232243",2378.09741210938 +"CR232244",3701.19702148438 +"CR232245",3164.34375 +"CR232246",2633.36059570312 +"CR232247",2589.1728515625 +"CR232248",2705.8232421875 +"CR232249",2573.6083984375 +"CR232250",2371.60424804688 +"CR232251",812.105163574219 +"CR232252",670.834106445312 +"CR232254",100.468406677246 +"CR232255",39.0627555847168 +"CR232256",12.9195594787598 +"CR233225",451.807373046875 +"CR233226",1497.94165039062 +"CR233227",815.475830078125 +"CR233229",1132.70397949219 +"CR233230",2327.62182617188 +"CR233231",1894.03601074219 +"CR233232",600.3212890625 +"CR233236",93.4322128295898 +"CR233237",64.6600646972656 +"CR233240",934.375549316406 +"CR233241",3040.46435546875 +"CR233242",4180.70458984375 +"CR233243",4602.578125 +"CR233244",4225.177734375 +"CR233245",4265.4267578125 +"CR233246",4090.42919921875 +"CR233247",3571.68237304688 +"CR233248",3392.01440429688 +"CR233249",2564.6982421875 +"CR233250",2249.11694335938 +"CR233251",2295.54858398438 +"CR233252",1880.69274902344 +"CR233253",984.16650390625 +"CR233254",659.370178222656 +"CR233255",343.208404541016 +"CR233256",45.4948654174805 +"CR234226",602.080139160156 +"CR234227",267.396423339844 +"CR234230",434.172515869141 +"CR234231",637.128479003906 +"CR234232",234.465896606445 +"CR234235",126.990875244141 +"CR234238",85.8135223388672 +"CR234240",529.359191894531 +"CR234241",809.768432617188 +"CR234242",3804.54174804688 +"CR234243",3680.26000976562 +"CR234244",3958.00659179688 +"CR234245",3220.25024414062 +"CR234246",3193.49682617188 +"CR234247",3668.1376953125 +"CR234248",3006.3107421875 +"CR234249",2633.87036132812 +"CR234250",2876.44604492188 +"CR234251",2860.65576171875 +"CR234252",2682.27514648438 +"CR234259",2.20393657684326 +"CR235232",45.1424903869629 +"CR235234",282.623016357422 +"CR235235",356.435272216797 +"CR235236",457.107177734375 +"CR235237",553.08203125 +"CR235238",1006.42114257812 +"CR235239",936.900024414062 +"CR235240",151.116912841797 +"CR235242",2701.69757952009 +"CR235243",2025.830078125 +"CR235244",2162.88208007812 +"CR235245",3754.37939453125 +"CR235246",3975.00830078125 +"CR235247",4023.73291015625 +"CR235248",4048.19848632812 +"CR235249",3223.12475585938 +"CR235250",2913.623046875 +"CR235251",2784.43676757812 +"CR235252",2759.42822265625 +"CR235258",39.5843048095703 +"CR235260",2.88824534416199 +"CR235262",3.5797438621521 +"CR236233",300.665313720703 +"CR236234",1018.86938476562 +"CR236235",968.990600585938 +"CR236236",1333.14428710938 +"CR236237",1060.49951171875 +"CR236238",2761.16967773438 +"CR236239",2147.55444335938 +"CR236240",785.169982910156 +"CR236241",668.490905761719 +"CR236242",1496.10327148438 +"CR236243",2066.65502929688 +"CR236244",2785.48657226562 +"CR236245",4055.71801757812 +"CR236246",4719.193359375 +"CR236247",4723.13623046875 +"CR236248",4458.71240234375 +"CR236249",3757.43603515625 +"CR236250",3053.64599609375 +"CR236251",2373.50659179688 +"CR236252",2434.02099609375 +"CR236253",2492.65600585938 +"CR236258",58.5163154602051 +"CR236260",1.20000004768372 +"CR237233",1236.82666015625 +"CR237234",2344.86767578125 +"CR237235",1745.79357910156 +"CR237236",1154.28820800781 +"CR237237",1080.1513671875 +"CR237238",1311.94860839844 +"CR237239",886.220336914062 +"CR237240",584.819213867188 +"CR237241",529.683410644531 +"CR237242",1774.77160644531 +"CR237243",3100.6826171875 +"CR237244",3744.93676757812 +"CR237245",4618.15185546875 +"CR237246",5121.755859375 +"CR237247",5072.4814453125 +"CR237248",4070.16650390625 +"CR237249",4438.31689453125 +"CR237250",3111.24243164062 +"CR237251",1932.35119628906 +"CR237252",2689.40966796875 +"CR237253",2616.126953125 +"CR237254",1663.89782714844 +"CR237255",848.58935546875 +"CR237258",334.615264892578 +"CR238234",2530.97998046875 +"CR238235",2882.76123046875 +"CR238236",623.652404785156 +"CR238237",638.389099121094 +"CR238238",638.821472167969 +"CR238239",595.933044433594 +"CR238240",333.836486816406 +"CR238241",457.476531982422 +"CR238242",1881.39086914062 +"CR238243",2820.57983398438 +"CR238244",3363.74340820312 +"CR238245",4423.708984375 +"CR238246",4934.0947265625 +"CR238247",5206.1025390625 +"CR238248",4498.63037109375 +"CR238249",5055.5234375 +"CR238250",3169.35009765625 +"CR238251",2711.10815429688 +"CR238252",2828.41430664062 +"CR238253",1965.81201171875 +"CR238254",1170.13330078125 +"CR238255",401.029937744141 +"CR239232",153.971115112305 +"CR239233",577.555603027344 +"CR239234",1617.56982421875 +"CR239235",1804.80627441406 +"CR239236",294.842315673828 +"CR239237",45.3207168579102 +"CR239238",66.0411758422852 +"CR239239",226.287231445312 +"CR239240",196.496841430664 +"CR239241",1029.09643554688 +"CR239242",2039.13317871094 +"CR239243",2179.57446289062 +"CR239244",2063.06201171875 +"CR239245",3313.47231445312 +"CR239246",3524.517578125 +"CR239247",5071.92138671875 +"CR239248",4666.4873046875 +"CR239249",4221.2080078125 +"CR239250",4039.37329101562 +"CR239251",3159.44506835938 +"CR239252",1876.15380859375 +"CR239253",1033.69567871094 +"CR239254",827.622680664062 +"CR239255",524.226806640625 +"CR239256",464.139617919922 +"CR240233",317.480438232422 +"CR240234",645.992797851562 +"CR240235",236.1337890625 +"CR240236",195.25862121582 +"CR240237",5.59088325500488 +"CR240238",44.9600105285645 +"CR240239",165.453582763672 +"CR240240",318.016418457031 +"CR240241",1357.27294921875 +"CR240242",1281.42053222656 +"CR240243",684.269287109375 +"CR240244",2718.86572265625 +"CR240245",3405.32934570312 +"CR240246",1006.44799804688 +"CR240247",2789.46606445312 +"CR240248",4767.42138671875 +"CR240249",4359.37060546875 +"CR240250",3787.85327148438 +"CR240251",1832.51525878906 +"CR240252",834.230834960938 +"CR240253",507.86376953125 +"CR240254",682.4443359375 +"CR240255",1046.52624511719 +"CR240256",1200.1806640625 +"CR240257",1298.54858398438 +"CR241232",31.1384944915771 +"CR241233",100.018966674805 +"CR241234",210.175720214844 +"CR241235",105.339859008789 +"CR241237",98.4239273071289 +"CR241238",548.829345703125 +"CR241239",657.655090332031 +"CR241240",175.313003540039 +"CR241241",256.03369140625 +"CR241242",966.351245117188 +"CR241243",1577.826171875 +"CR241244",1358.58337402344 +"CR241245",2386.52221679688 +"CR241246",2881.36694335938 +"CR241247",404.700866699219 +"CR241248",3160.43774414062 +"CR241249",3812.32958984375 +"CR241250",2229.16357421875 +"CR241251",1542.04174804688 +"CR241252",1367.85961914062 +"CR241253",515.418151855469 +"CR241254",840.059631347656 +"CR241255",1141.84167480469 +"CR241256",1157.86279296875 +"CR241257",1780.15869140625 +"CR242233",126.919128417969 +"CR242237",319.724029541016 +"CR242238",317.485260009766 +"CR242244",706.993774414062 +"CR242245",1172.60607910156 +"CR242246",2216.75634765625 +"CR242247",1497.81530761719 +"CR242248",729.153686523438 +"CR242249",2612.69311523438 +"CR242250",1265.27368164062 +"CR242251",1551.53918457031 +"CR242252",1255.95666503906 +"CR242253",618.862976074219 +"CR242254",602.5634765625 +"CR242255",1600.1650390625 +"CR242256",1530.24377441406 +"CR242257",1925.21118164062 +"CR242258",1698.39709472656 +"CR243232",28.4260425567627 +"CR243233",7.06616258621216 +"CR243234",92.8158569335938 +"CR243237",45.7404479980469 +"CR243245",2049.28564453125 +"CR243246",2028.28662109375 +"CR243248",64.9749755859375 +"CR243249",1023.13195800781 +"CR243250",678.706359863281 +"CR243251",1009.48754882812 +"CR243252",1008.02825927734 +"CR243253",333.931640625 +"CR243254",187.621704101562 +"CR243255",649.493103027344 +"CR243256",995.888122558594 +"CR243257",1115.90075683594 +"CR243258",1321.46105957031 +"CR244231",3.63072919845581 +"CR244233",1.20000004768372 +"CR244234",16.3856143951416 +"CR244237",15.8835773468018 +"CR244238",10.890585899353 +"CR244239",27.7999992370605 +"CR244245",1349.85437011719 +"CR244246",1712.07507324219 +"CR244249",957.5654296875 +"CR244250",401.790835571289 +"CR244251",200.760269165039 +"CR244252",356.366088867188 +"CR244253",44.8082275390625 +"CR244254",61.5291023254395 +"CR244255",77.8801956176758 +"CR244256",825.354858398438 +"CR244257",848.263671875 +"CR245231",1.8411111831665 +"CR245238",22.5169486999512 +"CR245239",106.180000305176 +"CR245240",95.7437744140625 +"CR245250",69.6325328826904 +"CR245251",71.9879302978516 +"CR245252",417.547271728516 +"CR245253",301.764434814453 +"CR245254",436.01416015625 +"CR245255",290.391723632812 +"CR245256",554.102478027344 +"CR245257",935.691589355469 +"CR246238",39.8641357421875 +"CR246251",172.284952163696 +"CR246252",254.591983795166 +"CR246253",257.54833984375 +"CR246254",335.326782226562 +"CR246255",444.581573486328 +"CR246256",244.496678670247 +"CR247253",161.393814086914 +"CR247254",136.650238037109 +"CR247255",181.127192179362 +"CR248235",856.519165039062 +"CR248236",109.688102722168 +"CR248237",44.7275276184082 +"CR249235",2721.71606445312 +"CR249236",1.69998979568481 +"CR249237",553.836669921875 +"CR250234",4106.9033203125 +"CR250235",1189.986328125 +"CR250236",166.417297363281 +"CR251232",4270.466796875 +"CR251234",1295.48742675781 +"CR251235",9.08649158477783 +"CR251236",2052.9970703125 +"CR252232",394.007751464844 +"CR252234",4.92215156555176 +"CR252235",760.611633300781 +"CR253232",2157.8388671875 +"CR253233",565.521240234375 +"CR253234",204.829467773438 +"CR155089",2.59999990463257 +"CR155090",2.52231860160828 +"CR155097",1.59993517398834 +"CR156083",89.9813613891602 +"CR156089",3.22327184677124 +"CR156090",2.15813779830933 +"CR156095",5.68578243255615 +"CR156097",1.89941585063934 +"CR157083",56.6881675720215 +"CR157084",90 +"CR157085",90 +"CR157087",966.5 +"CR157088",563.885803222656 +"CR157089",1.40037453174591 +"CR157091",8.8484468460083 +"CR157095",1.5574209690094 +"CR157096",4.23458242416382 +"CR157098",1.29999995231628 +"CR158083",57.6830024719238 +"CR158084",90 +"CR158085",90 +"CR158087",868.5 +"CR158088",734.067993164062 +"CR158089",2.44892072677612 +"CR158091",8.04464912414551 +"CR158095",3.2604341506958 +"CR158096",3.83363056182861 +"CR159083",70.0430755615234 +"CR159085",160.273162841797 +"CR159086",300 +"CR159087",920.5166015625 +"CR159088",283.298065185547 +"CR159089",4.39955139160156 +"CR159091",3.09999990463257 +"CR159095",3.42714405059814 +"CR159096",3.15921783447266 +"CR159097",8.99902057647705 +"CR159098",1.20000004768372 +"CR159099",1.10000002384186 +"CR160083",135.876739501953 +"CR160086",275.289367675781 +"CR160087",243.420471191406 +"CR160088",667.797180175781 +"CR160089",527.897583007812 +"CR160090",433.476348876953 +"CR160091",600.978149414062 +"CR160093",118.806137084961 +"CR160094",145.064590454102 +"CR160095",1.91675055027008 +"CR160096",2.64369702339172 +"CR160097",9.07360649108887 +"CR160098",2 +"CR160099",1.89999997615814 +"CR161083",240.283355712891 +"CR161084",579.800659179688 +"CR161085",993.071899414062 +"CR161086",506.918334960938 +"CR161087",228.703948974609 +"CR161088",689.531372070312 +"CR161089",672.295715332031 +"CR161090",743.263000488281 +"CR161091",390.7431640625 +"CR161092",228.427841186523 +"CR161093",232.108154296875 +"CR161094",220.983032226562 +"CR161096",7.62797021865845 +"CR161097",11.3999996185303 +"CR162083",388.057067871094 +"CR162084",592.948059082031 +"CR162085",937.47119140625 +"CR162086",1411.12036132812 +"CR162087",338.669738769531 +"CR162088",534.988830566406 +"CR162089",529.304809570312 +"CR162090",1089.5068359375 +"CR162091",526.628723144531 +"CR162092",276.732727050781 +"CR162093",732.66455078125 +"CR162094",382.215454101562 +"CR162095",56.6908950805664 +"CR162096",27.5320510864258 +"CR162097",4.40000009536743 +"CR162101",310 +"CR163083",279.799346923828 +"CR163084",616.115600585938 +"CR163085",870.866271972656 +"CR163086",2007.76281738281 +"CR163087",3428.52099609375 +"CR163088",1374.42883300781 +"CR163089",981.202819824219 +"CR163090",1355.33935546875 +"CR163091",1025.06420898438 +"CR163092",631.617126464844 +"CR163093",920.519897460938 +"CR163094",793.792236328125 +"CR163096",71.7851867675781 +"CR163101",460.799987792969 +"CR163102",136.309982299805 +"CR163107",175.716796875 +"CR164083",257.751190185547 +"CR164084",625.924011230469 +"CR164085",1020.35894775391 +"CR164086",1410.31384277344 +"CR164087",2385.22973632812 +"CR164088",1596.72351074219 +"CR164089",1328.14233398438 +"CR164090",1459.03149414062 +"CR164091",1502.78918457031 +"CR164092",952.972717285156 +"CR164093",1265.42663574219 +"CR164094",1384.44506835938 +"CR164095",960.376098632812 +"CR164102",456.651702880859 +"CR164103",260.455139160156 +"CR164106",155.993347167969 +"CR164107",129.252136230469 +"CR164108",186.575912475586 +"CR164112",39.4000015258789 +"CR165083",723.226318359375 +"CR165084",461.885131835938 +"CR165085",1388.85437011719 +"CR165086",1694.47473144531 +"CR165087",2011.84118652344 +"CR165088",1371.24182128906 +"CR165089",1161.36804199219 +"CR165090",1244.89306640625 +"CR165091",1334.67236328125 +"CR165092",1351.11682128906 +"CR165093",1484.03601074219 +"CR165094",1701.64489746094 +"CR165095",1238.24560546875 +"CR165102",561.944091796875 +"CR165103",478.507141113281 +"CR165106",399.348571777344 +"CR165107",273.296936035156 +"CR165108",163.050827026367 +"CR165117",2.4354612827301 +"CR166083",850.134582519531 +"CR166084",617.903991699219 +"CR166085",3733.40405273438 +"CR166086",3269.0712890625 +"CR166087",2992.68701171875 +"CR166088",1497.13647460938 +"CR166089",1047.82482910156 +"CR166090",1056.94921875 +"CR166091",1394.62512207031 +"CR166092",1494.19934082031 +"CR166093",1932.22265625 +"CR166094",1712.98889160156 +"CR166095",1679.87939453125 +"CR166101",1365.31567382812 +"CR166102",826.361633300781 +"CR166103",972.238037109375 +"CR166104",1177.07153320312 +"CR166105",956.434204101562 +"CR166106",874.856323242188 +"CR166107",901.2802734375 +"CR166108",178.278182983398 +"CR166109",114.015174865723 +"CR166111",59.1252784729004 +"CR166113",10.8016519546509 +"CR166114",4.78918313980103 +"CR166115",40.9136352539062 +"CR166116",14.9178438186646 +"CR166117",5.49405765533447 +"CR166118",45.6943817138672 +"CR166122",3.72614026069641 +"CR166123",57.5083236694336 +"CR167083",830.2001953125 +"CR167084",998.971130371094 +"CR167085",4025.34716796875 +"CR167086",2909.11499023438 +"CR167087",2249.87963867188 +"CR167088",2001.318359375 +"CR167089",1897.70458984375 +"CR167090",1423.08642578125 +"CR167091",1198.67211914062 +"CR167092",1295.21594238281 +"CR167093",1495.91015625 +"CR167094",2040.79150390625 +"CR167095",1850.65966796875 +"CR167099",1384.9228515625 +"CR167101",1450.17358398438 +"CR167102",128.054428100586 +"CR167103",724.094116210938 +"CR167104",415.545654296875 +"CR167105",265.508514404297 +"CR167106",299.988647460938 +"CR167107",144.350982666016 +"CR167108",109.048347473145 +"CR167109",59.8418197631836 +"CR167113",9.51197338104248 +"CR167114",9.40514373779297 +"CR167115",66.6556777954102 +"CR167116",57.3052978515625 +"CR167117",5.56277799606323 +"CR167118",34.3101501464844 +"CR167119",2.99732971191406 +"CR167120",4.00910139083862 +"CR167122",7.82928943634033 +"CR167123",201.373413085938 +"CR168083",587.175354003906 +"CR168084",1252.97387695312 +"CR168085",1569.00573730469 +"CR168086",2197.97485351562 +"CR168087",1488.61877441406 +"CR168088",1591.86657714844 +"CR168089",2093.77368164062 +"CR168090",1904.92614746094 +"CR168091",1928.15539550781 +"CR168092",1698.36767578125 +"CR168093",2367.41650390625 +"CR168094",2016.93615722656 +"CR168095",2088.43823242188 +"CR168096",1754.93115234375 +"CR168098",1693.52575683594 +"CR168099",1962.79370117188 +"CR168100",1142.15600585938 +"CR168101",1598.20031738281 +"CR168102",79.8900985717773 +"CR168103",541.835510253906 +"CR168104",642.694885253906 +"CR168105",288.649505615234 +"CR168106",165.820953369141 +"CR168107",39.8661575317383 +"CR168108",43.1454391479492 +"CR168109",41.0130081176758 +"CR168110",21.0364532470703 +"CR168111",7.72547101974487 +"CR168113",243.880813598633 +"CR168114",107.465751647949 +"CR168115",123.44108581543 +"CR168116",45.9260368347168 +"CR168118",2.32599973678589 +"CR168120",7.39009046554565 +"CR168121",8.27192211151123 +"CR168123",11.7984647750854 +"CR168124",50.838877360026 +"CR169083",608.687744140625 +"CR169084",1280.28784179688 +"CR169085",1106.56237792969 +"CR169086",902.504516601562 +"CR169087",436.363647460938 +"CR169088",849.760192871094 +"CR169089",1452.33947753906 +"CR169090",1627.5810546875 +"CR169091",2075 +"CR169092",1978.80236816406 +"CR169093",2251.06494140625 +"CR169094",2572.33837890625 +"CR169095",2587.8828125 +"CR169096",2046.69384765625 +"CR169097",1868.14807128906 +"CR169098",2052.70532226562 +"CR169099",1963.91442871094 +"CR169100",1893.06994628906 +"CR169101",1541.25952148438 +"CR169102",759.549987792969 +"CR169103",1180.7939453125 +"CR169104",1102.89904785156 +"CR169105",954.262084960938 +"CR169106",802.988342285156 +"CR169107",159.974716186523 +"CR169108",276.344177246094 +"CR169109",395.068420410156 +"CR169110",13.1000003814697 +"CR169113",410.122283935547 +"CR169114",463.021575927734 +"CR169115",150.699996948242 +"CR169122",21.3768405914307 +"CR169123",19.089225769043 +"CR170083",410.807037353516 +"CR170085",154.507598876953 +"CR170086",138.302093505859 +"CR170087",517.835144042969 +"CR170088",697.020690917969 +"CR170089",928.037292480469 +"CR170090",1392.53149414062 +"CR170091",1774.15991210938 +"CR170092",2123.6923828125 +"CR170093",2264.5888671875 +"CR170094",2120.69116210938 +"CR170095",2344.3857421875 +"CR170096",2012.93835449219 +"CR170097",2107.93383789062 +"CR170098",2149.74877929688 +"CR170099",1932.22351074219 +"CR170100",1888.568359375 +"CR170101",1752.04565429688 +"CR170102",905.235168457031 +"CR170103",1156.05749511719 +"CR170104",1199.521484375 +"CR170105",1392.13098144531 +"CR170106",1168.15466308594 +"CR170107",145.658721923828 +"CR170108",154.684387207031 +"CR170109",417.823516845703 +"CR170110",190.014266967773 +"CR170113",138.441482543945 +"CR170114",223.509216308594 +"CR170120",18.7736663818359 +"CR170122",15.3928689956665 +"CR171088",543.041015625 +"CR171089",477.808715820312 +"CR171090",1377.48693847656 +"CR171091",1730.12658691406 +"CR171092",1946.58166503906 +"CR171093",2338.68041992188 +"CR171094",1984.78283691406 +"CR171095",2408.3583984375 +"CR171096",2042.09118652344 +"CR171097",1941.56677246094 +"CR171098",1804.34680175781 +"CR171099",1625.99401855469 +"CR171100",1503.44555664062 +"CR171101",1165.57824707031 +"CR171102",908.701416015625 +"CR171103",783.651000976562 +"CR171104",905.087219238281 +"CR171105",750.524475097656 +"CR171106",1460.05346679688 +"CR171107",355.053680419922 +"CR171108",36.0369453430176 +"CR171109",38.8990440368652 +"CR171110",193.398696899414 +"CR171113",73.6301193237305 +"CR171114",61.0340538024902 +"CR171120",25.3064441680908 +"CR171121",80.5668411254883 +"CR172088",672.190673828125 +"CR172089",633.549133300781 +"CR172090",1147.08569335938 +"CR172091",1600.01086425781 +"CR172092",1886.99401855469 +"CR172093",2280.58325195312 +"CR172094",2504.42602539062 +"CR172095",2439.63598632812 +"CR172096",1934.31298828125 +"CR172097",1728.98498535156 +"CR172098",1480.19836425781 +"CR172099",995.537658691406 +"CR172100",1055.11608886719 +"CR172101",1233.322265625 +"CR172102",1140.39184570312 +"CR172103",1062.28137207031 +"CR172104",1499.62463378906 +"CR172105",504.376617431641 +"CR172107",71.402473449707 +"CR172109",223.68391418457 +"CR172110",89.0441131591797 +"CR172111",2.79999995231628 +"CR172113",273.758178710938 +"CR172114",11.3000001907349 +"CR172120",95.0999984741211 +"CR173089",1252.03466796875 +"CR173090",858.826904296875 +"CR173091",1043.28649902344 +"CR173092",1598.16271972656 +"CR173093",2286.82177734375 +"CR173094",1929.72424316406 +"CR173095",2177.92260742188 +"CR173096",1734.12573242188 +"CR173097",1380.99951171875 +"CR173098",1093.88012695312 +"CR173099",933.297424316406 +"CR173100",1408.29040527344 +"CR173101",1379.98999023438 +"CR173102",1946.31213378906 +"CR173103",1165.74609375 +"CR173104",613.924438476562 +"CR173105",264.785736083984 +"CR173106",337.710021972656 +"CR173107",75.8417892456055 +"CR173121",1.5 +"CR174088",217.592788696289 +"CR174089",492.987152099609 +"CR174090",681.082824707031 +"CR174091",469.394866943359 +"CR174092",1416.60034179688 +"CR174093",2004.19177246094 +"CR174094",1765.9951171875 +"CR174095",1812.3642578125 +"CR174096",1719.13598632812 +"CR174097",1727.67541503906 +"CR174098",1111.4755859375 +"CR174099",1010.25030517578 +"CR174100",1709.59338378906 +"CR174101",1833.07873535156 +"CR174102",1947.28381347656 +"CR174103",1376.90051269531 +"CR174104",286.870574951172 +"CR174105",66.7453079223633 +"CR174106",64.3085861206055 +"CR174115",2.20000004768372 +"CR174118",13.1000003814697 +"CR174119",72.0999984741211 +"CR175088",74.0999984741211 +"CR175089",259.664459228516 +"CR175090",672.197998046875 +"CR175091",659.58251953125 +"CR175092",1295.8544921875 +"CR175093",1829.29858398438 +"CR175094",2011.25756835938 +"CR175095",2079.02685546875 +"CR175096",2074.73461914062 +"CR175097",1956.23876953125 +"CR175098",1495.30322265625 +"CR175099",847.611267089844 +"CR175100",877.081420898438 +"CR175101",1465.77709960938 +"CR175102",1783.42602539062 +"CR175103",1093.19458007812 +"CR175104",140.55696105957 +"CR175110",145.812286376953 +"CR175111",136.556503295898 +"CR175117",176.625381469727 +"CR175118",336.9208984375 +"CR175119",29.9058418273926 +"CR176088",68.7834930419922 +"CR176089",86.4593200683594 +"CR176090",485.430114746094 +"CR176091",596.702514648438 +"CR176092",699.904907226562 +"CR176093",1216.01000976562 +"CR176094",1657.82189941406 +"CR176095",1623.88110351562 +"CR176096",1814.509765625 +"CR176097",1695.48937988281 +"CR176098",1569.80627441406 +"CR176099",1007.05920410156 +"CR176100",1247.05151367188 +"CR176101",1734.234375 +"CR176102",1329.84643554688 +"CR176103",531.8798828125 +"CR176104",172.764266967773 +"CR176110",165.511276245117 +"CR176111",51.4722099304199 +"CR176112",23.7931823730469 +"CR176117",9.58837699890137 +"CR176118",590.402526855469 +"CR176119",919.272827148438 +"CR176120",1221.90234375 +"CR177088",36.5929718017578 +"CR177089",132.750457763672 +"CR177090",524.205993652344 +"CR177091",553.487426757812 +"CR177092",364.425537109375 +"CR177093",817.445495605469 +"CR177094",869.115478515625 +"CR177095",1332.51416015625 +"CR177096",1360.98901367188 +"CR177097",1220.43273925781 +"CR177098",1643.66345214844 +"CR177099",1120.01586914062 +"CR177100",1329.35339355469 +"CR177101",1464.90417480469 +"CR177102",1695.69775390625 +"CR177103",495.069396972656 +"CR177104",132.200942993164 +"CR177109",140.300964355469 +"CR177110",159.927444458008 +"CR177114",7.69999980926514 +"CR177118",777.434936523438 +"CR177119",1157.88073730469 +"CR177120",626.531982421875 +"CR178089",68.4698257446289 +"CR178090",230.536178588867 +"CR178091",371.831726074219 +"CR178092",467.182312011719 +"CR178093",352.247741699219 +"CR178094",468.91650390625 +"CR178095",683.330383300781 +"CR178096",1047.35913085938 +"CR178097",1749.98254394531 +"CR178098",1580.47106933594 +"CR178099",1300.04772949219 +"CR178100",1591.74230957031 +"CR178101",1342.89990234375 +"CR178102",1315.80969238281 +"CR178103",791.348327636719 +"CR178104",231.05615234375 +"CR178107",53.0014228820801 +"CR178109",207.199996948242 +"CR178117",116.345634460449 +"CR178119",926.079895019531 +"CR178120",227.573028564453 +"CR178121",54.2792701721191 +"CR179089",80 +"CR179090",130.824752807617 +"CR179091",252.489242553711 +"CR179092",191.671493530273 +"CR179093",366.152008056641 +"CR179094",396.640930175781 +"CR179095",470.88818359375 +"CR179096",625.684326171875 +"CR179097",1704.47692871094 +"CR179098",1438.07592773438 +"CR179099",1784.86889648438 +"CR179100",1368.50048828125 +"CR179101",1213.15478515625 +"CR179102",911.161865234375 +"CR179103",730.587585449219 +"CR179104",165.862762451172 +"CR179106",34.4558181762695 +"CR179107",182.608383178711 +"CR179116",176.757583618164 +"CR180090",140.989334106445 +"CR180091",300 +"CR180092",212.481842041016 +"CR180093",347.036499023438 +"CR180094",307.687408447266 +"CR180095",651.143371582031 +"CR180096",534.101989746094 +"CR180097",1461.91381835938 +"CR180098",1784.31945800781 +"CR180099",1558.85522460938 +"CR180100",1197.66955566406 +"CR180101",1645.17602539062 +"CR180102",1665.01220703125 +"CR180103",857.75537109375 +"CR180104",1421.93115234375 +"CR180105",175.241836547852 +"CR180106",29.8903980255127 +"CR180107",520.908935546875 +"CR180110",348.198364257812 +"CR180115",381.794982910156 +"CR180116",111.154411315918 +"CR181090",36.3629531860352 +"CR181091",970.817565917969 +"CR181092",1000.27532958984 +"CR181093",353.512573242188 +"CR181094",405.936370849609 +"CR181095",624.416137695312 +"CR181096",1066.08129882812 +"CR181097",1516.43322753906 +"CR181098",1840.34997558594 +"CR181099",1129.90417480469 +"CR181100",1794.91796875 +"CR181101",1432.064453125 +"CR181102",1730.88916015625 +"CR181103",2098.69848632812 +"CR181104",2062.05224609375 +"CR181105",1235.32397460938 +"CR181106",1140.53247070312 +"CR181108",3142.1748046875 +"CR181109",2239.01586914062 +"CR181110",932.256958007812 +"CR181111",329.834289550781 +"CR181112",137.024108886719 +"CR181113",125.95832824707 +"CR181114",59.71826171875 +"CR181115",86.459098815918 +"CR181116",109.318130493164 +"CR182089",40 +"CR182090",40 +"CR182091",189.8984375 +"CR182092",597.410766601562 +"CR182093",442.379486083984 +"CR182094",547.197326660156 +"CR182095",908.167541503906 +"CR182096",931.565551757812 +"CR182097",1251.03381347656 +"CR182098",1404.54736328125 +"CR182099",1810.26135253906 +"CR182100",2207.16943359375 +"CR182101",1894.53039550781 +"CR182102",2157.6123046875 +"CR182103",2193.49877929688 +"CR182104",2063.01123046875 +"CR182105",1165.07043457031 +"CR182106",1204.98815917969 +"CR182107",2139.11010742188 +"CR182108",2403.73706054688 +"CR182109",1837.55029296875 +"CR182110",641.961547851562 +"CR182111",340.244445800781 +"CR182112",350.014556884766 +"CR182113",236.896408081055 +"CR183090",178.730667114258 +"CR183091",295.243347167969 +"CR183092",710.3193359375 +"CR183093",697.694458007812 +"CR183094",787.337646484375 +"CR183095",1052.67211914062 +"CR183096",1222.45373535156 +"CR183097",1546.7314453125 +"CR183098",1761.40979003906 +"CR183099",2391.087890625 +"CR183100",2260.71997070312 +"CR183101",1935.421875 +"CR183102",2082.33862304688 +"CR183103",1857.67565917969 +"CR183104",1626.78735351562 +"CR183105",961.444152832031 +"CR183106",1025.90612792969 +"CR183107",1655.224609375 +"CR183108",1654.27185058594 +"CR183109",845.46630859375 +"CR183110",345.772216796875 +"CR183111",445.341400146484 +"CR183112",826.676513671875 +"CR184089",49.9914665222168 +"CR184090",95.5637435913086 +"CR184091",244.864456176758 +"CR184092",738.674438476562 +"CR184093",519.112854003906 +"CR184094",414.611450195312 +"CR184095",886.283203125 +"CR184096",857.430236816406 +"CR184097",968.326599121094 +"CR184098",1640.03112792969 +"CR184099",2519.72485351562 +"CR184100",2275.94384765625 +"CR184101",2250.27685546875 +"CR184102",2107.90649414062 +"CR184103",2500.57006835938 +"CR184104",2349.669921875 +"CR184105",1732.08386230469 +"CR184106",597.089172363281 +"CR184107",1051.85487060547 +"CR184108",816.520202636719 +"CR184109",343.364440917969 +"CR184110",431.871368408203 +"CR184111",123.445663452148 +"CR184112",113.64714050293 +"CR184113",8.40706157684326 +"CR184114",56.6447296142578 +"CR184115",22.9756126403809 +"CR184116",5.12067365646362 +"CR185090",50.0269546508789 +"CR185091",235.046539306641 +"CR185092",325.329559326172 +"CR185093",396.351318359375 +"CR185094",583.671821594238 +"CR185095",804.619247436523 +"CR185096",790.235681152344 +"CR185097",658.248107910156 +"CR185098",1416.02368164062 +"CR185099",2264.9970703125 +"CR185100",2284.21508789062 +"CR185101",2254.54223632812 +"CR185102",1938.95959472656 +"CR185103",2333.32275390625 +"CR185104",1912.85754394531 +"CR185105",1790.35754394531 +"CR185106",751.059875488281 +"CR185107",626.638549804688 +"CR185108",88.7549362182617 +"CR185109",88.4778060913086 +"CR185110",57.903980255127 +"CR185111",291.466033935547 +"CR185112",35.6325950622559 +"CR185113",45.685905456543 +"CR185114",64.8405303955078 +"CR185115",66.9234619140625 +"CR185118",81.5090560913086 +"CR185119",157.921463012695 +"CR185120",417.836181640625 +"CR186090",73.9919662475586 +"CR186091",157.417922973633 +"CR186097",914.545837402344 +"CR186098",1551.87670898438 +"CR186099",2278.54370117188 +"CR186100",1873.29504394531 +"CR186101",1428.564453125 +"CR186102",1783.0302734375 +"CR186103",1662.71130371094 +"CR186104",1522.16979980469 +"CR186105",1787.00842285156 +"CR186106",775.411010742188 +"CR186107",1187.99877929688 +"CR186108",163.462753295898 +"CR186109",134.555557250977 +"CR186110",290.123565673828 +"CR186111",512.597351074219 +"CR186112",212.030899047852 +"CR186113",25.1460342407227 +"CR186114",14.8999996185303 +"CR186116",99.2472076416016 +"CR186118",87.1292572021484 +"CR186119",41.4754409790039 +"CR186120",197.703842163086 +"CR187097",427.715515136719 +"CR187098",1872.39501953125 +"CR187099",1901.068359375 +"CR187100",2008.69445800781 +"CR187101",1842.74987792969 +"CR187102",656.964111328125 +"CR187103",682.550598144531 +"CR187104",821.891967773438 +"CR187105",622.589599609375 +"CR187106",428.895172119141 +"CR187107",1316.87048339844 +"CR187108",263.339691162109 +"CR187109",182.956985473633 +"CR187110",232.752166748047 +"CR187111",664.32275390625 +"CR187112",120.531204223633 +"CR187113",26.0070571899414 +"CR187116",113.267097473145 +"CR187119",11.1219253540039 +"CR188092",52.0589141845703 +"CR188093",61.011287689209 +"CR188096",499.986145019531 +"CR188098",1377.99621582031 +"CR188099",1966.11572265625 +"CR188100",2213.10620117188 +"CR188101",918.398681640625 +"CR188102",604.727172851562 +"CR188103",1373.70068359375 +"CR188104",749.000793457031 +"CR188105",478.860229492188 +"CR188106",292.720336914062 +"CR188107",327.573883056641 +"CR188108",132.530303955078 +"CR188109",332.983001708984 +"CR188110",463.399261474609 +"CR188111",321.536682128906 +"CR188112",103.988731384277 +"CR188113",98.602783203125 +"CR188119",19.6899852752686 +"CR189095",495.281311035156 +"CR189098",1794.37963867188 +"CR189099",2337.8603515625 +"CR189100",2215.39208984375 +"CR189101",2286.86572265625 +"CR189102",2047.08276367188 +"CR189103",1347.57104492188 +"CR189104",527.526489257812 +"CR189105",469.128082275391 +"CR189106",299.648681640625 +"CR189107",122.49568939209 +"CR189108",51.797779083252 +"CR189109",264.187591552734 +"CR189110",275.597015380859 +"CR189111",281.298919677734 +"CR189112",245.031509399414 +"CR189113",92.7322540283203 +"CR189119",22.5606861114502 +"CR190093",476.964141845703 +"CR190094",713.57666015625 +"CR190095",780.050048828125 +"CR190096",863.368896484375 +"CR190097",1155.05126953125 +"CR190098",1552.15649414062 +"CR190099",2764.01025390625 +"CR190100",2282.35791015625 +"CR190101",1912.27917480469 +"CR190102",1380.37390136719 +"CR190103",688.701049804688 +"CR190104",570.963317871094 +"CR190105",334.368377685547 +"CR190106",229.975860595703 +"CR190107",56.6401748657227 +"CR190108",23.3479995727539 +"CR190109",34.822940826416 +"CR190110",90.5826797485352 +"CR190111",45.3726615905762 +"CR190112",82.7936477661133 +"CR190119",56.1678199768066 +"CR191093",847.516418457031 +"CR191094",1144.7578125 +"CR191095",1178.45251464844 +"CR191096",1014.51696777344 +"CR191097",1200.28479003906 +"CR191098",2451.34423828125 +"CR191099",2375.77709960938 +"CR191100",2224.546875 +"CR191101",1910.97888183594 +"CR191102",600.336547851562 +"CR191103",299.914703369141 +"CR191105",135.053176879883 +"CR191106",178.061111450195 +"CR191107",116.353073120117 +"CR191109",105.522690582275 +"CR191110",118.019844055176 +"CR191111",164.698272705078 +"CR191112",81.1568603515625 +"CR191116",49.0953750610352 +"CR191119",17.7455711364746 +"CR192092",25.7153511047363 +"CR192093",632.053527832031 +"CR192094",1303.27648925781 +"CR192095",1020.01495361328 +"CR192096",659.706726074219 +"CR192097",1700.49951171875 +"CR192098",2430.7099609375 +"CR192099",2571.01416015625 +"CR192100",2102.63061523438 +"CR192101",1109.92211914062 +"CR192102",853.381591796875 +"CR192103",882.792724609375 +"CR192104",263.804351806641 +"CR192105",63.6336784362793 +"CR192106",56.2621994018555 +"CR192107",1.5 +"CR192109",162.43586730957 +"CR192110",53.2432861328125 +"CR192115",193.220825195312 +"CR192116",115.828994750977 +"CR192119",7.0899486541748 +"CR193092",32.7732353210449 +"CR193093",1054.68179931641 +"CR193094",1134.61950683594 +"CR193095",710.623168945312 +"CR193096",836.102478027344 +"CR193097",1370.27270507812 +"CR193098",2274.61303710938 +"CR193099",2337.49658203125 +"CR193100",1941.8017578125 +"CR193101",2265.92163085938 +"CR193102",2425.33447265625 +"CR193103",866.718139648438 +"CR193104",164.784744262695 +"CR193105",103.741554260254 +"CR193109",50.2451667785645 +"CR193114",28.348123550415 +"CR193115",52.0233383178711 +"CR193116",113.699996948242 +"CR193119",10.1999998092651 +"CR194093",1733.10522460938 +"CR194094",1094.02673339844 +"CR194095",634.711669921875 +"CR194096",1069.07348632812 +"CR194097",1660.30491420201 +"CR194098",2318.07934570312 +"CR194099",2303.27416992188 +"CR194100",2100.98364257812 +"CR194101",2327.7236328125 +"CR194102",1034.66748046875 +"CR194103",251.225234985352 +"CR194104",27.891508102417 +"CR194109",16.8455982208252 +"CR194111",46.1496391296387 +"CR194114",33.388484954834 +"CR194115",240.701934814453 +"CR194119",6.69999980926514 +"CR195093",954.353881835938 +"CR195094",1567.36389160156 +"CR195095",843.704284667969 +"CR195098",1964.13305664062 +"CR195099",1733.99499511719 +"CR195100",1366.15625 +"CR195101",1334.92846679688 +"CR195102",405.915649414062 +"CR195103",170.817031860352 +"CR195104",24.9169502258301 +"CR195108",38.4057426452637 +"CR195111",50.1430320739746 +"CR195112",39.6940307617188 +"CR195114",51.5060348510742 +"CR195115",217.826248168945 +"CR195116",206.249603271484 +"CR196094",1701.25687953404 +"CR196095",2254 +"CR196098",1301.984375 +"CR196099",852.080932617188 +"CR196100",694.8173828125 +"CR196101",470.851196289062 +"CR196102",43.7837600708008 +"CR196103",94.8886260986328 +"CR196104",56.8327713012695 +"CR196111",30.4364700317383 +"CR196112",44.0999984741211 +"CR196115",90.3981781005859 +"CR196116",144.380493164062 +"CR196117",144.642349243164 +"CR196118",34.5909729003906 +"CR197098",450.949798583984 +"CR197099",638.89501953125 +"CR197100",308.058044433594 +"CR197101",113.105369567871 +"CR197102",74.9365921020508 +"CR197103",149.656814575195 +"CR197109",136.110824584961 +"CR197110",241.072296142578 +"CR197111",34.5999984741211 +"CR197112",104.956634521484 +"CR197113",72.0636291503906 +"CR197114",49.5837669372559 +"CR197115",346.521575927734 +"CR197116",242.163833618164 +"CR197118",15.3352537155151 +"CR198098",333.191223144531 +"CR198099",384.148162841797 +"CR198100",149.190582275391 +"CR198101",63.6737365722656 +"CR198102",72.1857986450195 +"CR198109",200.334579467773 +"CR198110",187.923004150391 +"CR198111",131.257415771484 +"CR198112",37.3310928344727 +"CR198113",96.3323745727539 +"CR198114",586.374877929688 +"CR198115",277.012329101562 +"CR198116",64.4634628295898 +"CR199097",965.089965820312 +"CR199098",691.137268066406 +"CR199099",490.442565917969 +"CR199100",50.7056655883789 +"CR199107",84.4970092773438 +"CR199108",376.534484863281 +"CR199109",295.325653076172 +"CR199110",410.381317138672 +"CR199111",359.618255615234 +"CR199112",106.649063110352 +"CR199113",170.678756713867 +"CR199114",215.178802490234 +"CR199115",71.5838623046875 +"CR200096",191.514101664225 +"CR200097",227.5322265625 +"CR200098",256.038024902344 +"CR200099",103.189910888672 +"CR200100",119.714202880859 +"CR200106",20.7901058197021 +"CR200107",35.6573219299316 +"CR200108",107.611015319824 +"CR200109",190.960815429688 +"CR200110",528.627990722656 +"CR200111",300.491607666016 +"CR200112",727.392456054688 +"CR200113",548.363952636719 +"CR200114",70.5830230712891 +"CR200115",30.7233352661133 +"CR201096",133.535726928711 +"CR201097",99.8137359619141 +"CR201098",45.4063301086426 +"CR201099",80.1422958374023 +"CR201100",73.0300521850586 +"CR201101",41.7854194641113 +"CR201106",27.2733535766602 +"CR201107",45.2961387634277 +"CR201108",148.650650024414 +"CR201109",168.766357421875 +"CR201110",41.3074340820312 +"CR201111",398.99951171875 +"CR201112",606.46484375 +"CR201113",628.181213378906 +"CR201114",93.989143371582 +"CR201115",13.1146806081136 +"CR202094",1666.73701171875 +"CR202095",1007.28431828817 +"CR202098",39.8739395141602 +"CR202099",114.003448486328 +"CR202100",73.2214584350586 +"CR202101",48.1695747375488 +"CR202103",5.37738084793091 +"CR202104",7.75091600418091 +"CR202105",32.5430068969727 +"CR202106",45.831111907959 +"CR202107",56.8421821594238 +"CR202108",176.607879638672 +"CR202109",78.310173034668 +"CR202110",178.017532348633 +"CR202111",996.202453613281 +"CR202112",743.004150390625 +"CR202113",281.648651123047 +"CR203094",303.8203125 +"CR203095",98.5496444702148 +"CR203099",64.5528945922852 +"CR203100",46.2022476196289 +"CR203101",34.8503608703613 +"CR203103",45.1280288696289 +"CR203104",20.4606151580811 +"CR203105",6.9148211479187 +"CR203106",71.5107192993164 +"CR203107",110.761184692383 +"CR203108",191.575866699219 +"CR203109",295.736694335938 +"CR203110",567.323669433594 +"CR203111",461.873809814453 +"CR203112",664.012817382812 +"CR203113",214.640106201172 +"CR204094",430.00537109375 +"CR204095",297.958862304688 +"CR204098",7.02562570571899 +"CR204099",37.0198745727539 +"CR204100",54.4711494445801 +"CR204101",37.6108856201172 +"CR204102",60.8355598449707 +"CR204103",58.7042808532715 +"CR204104",128.885696411133 +"CR204105",26.570384979248 +"CR204106",83.8175201416016 +"CR204107",180.240844726562 +"CR204108",330.2822265625 +"CR204109",858.899780273438 +"CR204110",931.687316894531 +"CR204111",439.203552246094 +"CR204112",202.551620483398 +"CR205094",196.271179199219 +"CR205095",441.509307861328 +"CR205098",37.8721542358398 +"CR205099",140.447845458984 +"CR205100",109.833396911621 +"CR205101",265.588653564453 +"CR205102",301.190795898438 +"CR205103",147.565277099609 +"CR205104",126.864456176758 +"CR205105",161.115295410156 +"CR205106",223.747955322266 +"CR205107",302.659820556641 +"CR205108",468.545013427734 +"CR205109",914.519714355469 +"CR205110",946.376647949219 +"CR205111",461.675170898438 +"CR205112",151.050542195638 +"CR206094",349.413500976562 +"CR206095",364.003204345703 +"CR206098",23.0143623352051 +"CR206099",204.777099609375 +"CR206100",271.2548828125 +"CR206101",475.339416503906 +"CR206102",428.571228027344 +"CR206103",53.136531829834 +"CR206104",322.384216308594 +"CR206105",499.667633056641 +"CR206106",368.83173828125 +"CR206107",538.733215332031 +"CR206108",505.977142333984 +"CR206109",756.0361328125 +"CR206110",577.150817871094 +"CR206111",167.286026000977 +"CR207094",424.825775146484 +"CR207095",556.258972167969 +"CR207096",142.197540283203 +"CR207098",499.586669921875 +"CR207099",338.204742431641 +"CR207100",281.501159667969 +"CR207101",474.721893310547 +"CR207102",148.18489074707 +"CR207103",162.93212890625 +"CR207104",479.263488769531 +"CR207105",491.119000244141 +"CR207106",477.038269042969 +"CR207107",633.935485839844 +"CR207108",587.613464355469 +"CR207109",841.723510742188 +"CR207110",470.904968261719 +"CR207111",434.748687744141 +"CR208092",421.167209201389 +"CR208093",82.8175048828125 +"CR208094",126.641441345215 +"CR208095",143.974990844727 +"CR208096",33.5159378051758 +"CR208098",127.476409912109 +"CR208099",129.098129272461 +"CR208100",416.349182128906 +"CR208101",362.996459960938 +"CR208102",447.759666442871 +"CR208103",1082.47560628255 +"CR208104",344.597310384115 +"CR208105",584.938339233398 +"CR208106",508.311367034912 +"CR208107",893.447021484375 +"CR208108",1641.32800292969 +"CR208109",674.48681640625 +"CR208110",384.82216389974 +"CR208111",546.928558349609 +"CR209092",88.9275512695312 +"CR209093",63.0999984741211 +"CR209094",143.933792114258 +"CR209095",33.8627624511719 +"CR209098",23.0268535614014 +"CR209099",281.827758789062 +"CR209100",217.666702270508 +"CR209101",284.991760253906 +"CR209102",1159.73413085938 +"CR209103",2284.99682617188 +"CR209104",821.002807617188 +"CR209105",907.569946289062 +"CR209106",769.220123291016 +"CR209107",754.8212890625 +"CR209108",758.641276041667 +"CR210094",131.326934814453 +"CR210095",14.9572172164917 +"CR210098",10.5341215133667 +"CR210099",237.789169311523 +"CR210100",213.606796264648 +"CR210101",757.072631835938 +"CR210102",321.700012207031 +"CR210103",16.8999996185303 +"CR210104",875.353454589844 +"CR210105",558.019327799479 +"CR211091",15.6859731674194 +"CR211095",34.7033843994141 +"CR211098",58.4474945068359 +"CR211099",105.734931945801 +"CR211100",248.699432373047 +"CR211101",241.176513671875 +"CR212095",76.4141006469727 +"CR212099",13.5 +"CR212100",9.33709049224854 +"CR212101",18.7145347595215 +"CR213096",63.3865089416504 diff --git a/code/tiff_to_nc.py b/code/tiff_to_nc.py new file mode 100644 index 0000000..914f0aa --- /dev/null +++ b/code/tiff_to_nc.py @@ -0,0 +1,592 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Jun 22 16:05:19 2021 + +@author: morenodu +""" +import os +os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data') +import xarray as xr +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd +import geopandas as gpd +import cartopy +import cartopy.crs as ccrs +import cartopy.io.shapereader as shpreader +from scipy import signal +import seaborn as sns +import matplotlib as mpl +from osgeo import gdal +import xarray as xr +import matplotlib.pyplot as plt +import cartopy.crs as ccrs + +#Change the following variables to the file you want to convert (inputfile) and +#what you want to name your output file (outputfile). +inputfile = 'spam2010v1r0_global_harvested-area_soyb_r.tif' +outputfile = 'spam2010v1r0_global_harvested-area_soyb_r.nc' +#Do not change this line, the following command will convert the geoTIFF to a netCDF +ds = gdal.Translate(outputfile, inputfile, format='NetCDF') + +DS_y_obs = xr.open_dataset("spam2010v1r0_global_harvested-area_soyb_r.nc", decode_times=True) + +def plot_2d_map(dataarray_2d): + # Plot 2D map of DataArray, remember to average along time or select one temporal interval + plt.figure(figsize=(12,5)) #plot clusters + ax=plt.axes(projection=ccrs.Mercator()) + dataarray_2d.plot(x='lon', y='lat',transform=ccrs.PlateCarree(), robust=True) + # ax.add_geometries(us1_shapes, ccrs.PlateCarree(),edgecolor='black', facecolor=(0,1,0,0.0)) + ax.set_extent([-80.73,-34,-35,6], ccrs.PlateCarree()) + plt.show() + +plot_2d_map(DS_y_obs.Band1) + + + +DS_y_epic = xr.open_dataset("epic-iiasa_gswp3-w5e5_obsclim_2015soc_default_yield-soy-noirr_global_annual_1901_2016.nc", decode_times=False) +DS_y_epic_mean = DS_y_epic.mean('time') +DS_y_epic_mean.to_netcdf("epic_grid_05x05.nc") + +#%% TRANSFORM HARVEST AREA + +# Load Geotiff +os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data') + +# US +rio = xr.open_rasterio("soy_harvest_area_us_1980_2016_05x05_density_02.tif") +rio = rio.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +rio['time'] = list(map(int, rio.attrs['descriptions'])) +rio.name = 'harvest_area' +# rio = rio/1000 + +ds_rio = rio.to_dataset() +ds_rio['time'] = pd.date_range(start='1980', periods=ds_rio.sizes['time'], freq='YS').year +# ds_rio = ds_rio.mean('time') +ds_rio['time'].attrs = {'units':'years since 1980-01-01'} +ds_rio['lat'].attrs = {'units':'degrees_north'} +ds_rio['lon'].attrs = {'units':'degrees_east'} +# ds_rio = ds_rio.where(ds_rio['harvest_area'] > -1000) + +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio['harvest_area'].sum(['lat','lon']).plot(label = 'US') +ds_rio['harvest_area'].sum(['lat','lon']).sel(time = 2010) + +ds_rio.to_netcdf("soy_harvest_area_US_all_1980_2020_05x05_density.nc") + +plot_2d_am_map(ds_rio['harvest_area'].mean('time')) + +df_rio = rio.to_dataframe().dropna() +df_rio_mean = df_rio.groupby('time').mean(...) + +#### +# ARG +rio_2 = xr.open_rasterio("soy_harvest_area_arg_1978_2000_05x05_density_02.tif") +rio_2 = rio_2.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +rio['time'] = list(map(int, rio.attrs['descriptions'])) +rio_2.name = 'harvest_area' +# rio = rio/1000 + +ds_rio = rio_2.to_dataset() +ds_rio['time'] = pd.date_range(start='1978', periods=ds_rio.sizes['time'], freq='YS').year +ds_rio['time'].attrs = {'units':'years since 1978-01-01'} +# ds_rio = ds_rio.mean('time') +ds_rio['lat'].attrs = {'units':'degrees_north'} +ds_rio['lon'].attrs = {'units':'degrees_east'} + +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio.to_netcdf("soy_harvest_area_arg_1978_2019_05x05_density.nc") + +plot_2d_am_map(ds_rio['harvest_area'].sel(time=1990)) + +df_rio = rio.to_dataframe().dropna() +df_rio_mean = df_rio.groupby('time').mean(...) + + +# Brazil +rio_2 = xr.open_rasterio("soy_harvest_area_br_1980_2016_05x05_density_03.tif") +rio_2 = rio_2.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +rio_2['time'] = list(map(int, rio_2.attrs['descriptions'])) +rio_2.name = 'harvest_area' +# rio = rio/1000 + +ds_rio = rio_2.to_dataset() +ds_rio['time'] = pd.date_range(start='1980', periods=ds_rio.sizes['time'], freq='YS').year +# ds_rio = ds_rio.mean('time') +ds_rio['lat'].attrs = {'units':'degrees_north'} +ds_rio['lon'].attrs = {'units':'degrees_east'} +ds_rio['time'].attrs = {'units':'years since 1980-01-01'} + +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio.to_netcdf("soy_harvest_area_br_1980_2016_05x05_density.nc") + +plot_2d_am_map(ds_rio['harvest_area'].sel(time=2016)) + +df_rio = rio.to_dataframe().dropna() +df_rio_mean = df_rio.groupby('time').mean(...) + +# SPAM + +rio_2 = xr.open_rasterio("soy_harvest_spam_agg_resamp.tif") +rio_2 = rio_2.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) + +rio_2.name = 'harvest_area' +rio_2 = rio_2.mean('time') +ds_rio = rio_2.to_dataset() +plot_2d_am_map(ds_rio['harvest_area'].where(ds_rio['harvest_area']>0)) + +ds_rio.to_netcdf("soy_harvest_spam_agg_resamp.nc") + +# SPAM 2 +rio_spam = xr.open_rasterio("soy_harvest_spam_native.tif") +rio_spam = rio_spam.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) + +rio_spam.name = 'harvest_area' +rio_spam = rio_spam.mean('time') +ds_rio_spam = rio_spam.to_dataset() +ds_rio_spam['lat'].attrs = {'units':'degrees_north'} +ds_rio_spam['lon'].attrs = {'units':'degrees_east'} +plot_2d_am_map(ds_rio_spam['harvest_area'].where(ds_rio_spam['harvest_area']>10)) + +ds_rio_spam.to_netcdf("soy_harvest_spam_native.nc") + +#%% TRANSFORM YIELD + +# Load Geotiff +os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data') + +rio = xr.open_rasterio("soy_yields_US_all_1975_2020_1prc_002.tif") +rio = rio.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +# rio['time'] = list(map(int, rio.attrs['descriptions'])) +rio.name = 'Yield' +# rio = rio/1000 + +ds_rio = rio.to_dataset() +# ds_rio = ds_rio.where(ds_rio['harvest_area'] > -1000) +ds_rio['time'] = pd.date_range(start='1975', periods=ds_rio.sizes['time'], freq='YS').year + +ds_rio['time'].attrs = {'units':'years since 1975-01-01'} +ds_rio['lat'].attrs = {'units':'degrees_north'} +ds_rio['lon'].attrs = {'units':'degrees_east'} +ds_rio = ds_rio.sortby('lat') +ds_rio = ds_rio.sortby('lon') +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio.to_netcdf("soy_yields_US_all_1975_2020_1prc_002.nc") + +plot_2d_am_map(ds_rio['Yield'].sel(time = 1987)) +plot_2d_am_map(ds_rio['Yield'].sel(time = 1990)) +plot_2d_am_map(ds_rio['Yield'].mean('time')) + + +# ARGE +rio_arg = xr.open_rasterio("soy_yield_arg_1978_2019_005.tif") +rio_arg = rio_arg.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +rio_arg['time'] = list(map(int, rio_arg.attrs['descriptions'])) +rio_arg.name = 'Yield' + +ds_rio_arg = rio_arg.to_dataset() +# ds_rio = ds_rio.where(ds_rio['harvest_area'] > -1000) +ds_rio_arg = ds_rio_arg.sel(lat = slice(-15,-55), lon = slice(-90,-40)) # SLICE IT TO BECOME FASTER TO CHANGE RESOLUTION LATER ON CDO +ds_rio_arg['Yield'] = ds_rio_arg['Yield']/1000 + +ds_rio_arg['lat'].attrs = {'units':'degrees_north'} +ds_rio_arg['lon'].attrs = {'units':'degrees_east'} +ds_rio_arg['time'].attrs = {'units':f'years since {rio_arg.time[0]-1}-01-01'} +ds_rio_arg = ds_rio_arg.sortby('lat') +ds_rio_arg = ds_rio_arg.sortby('lon') +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio_arg.to_netcdf("soy_yield_arg_1974_2019_005.nc") + +plot_2d_am_map(ds_rio_arg['Yield'].mean('time')) + + + +# BRASIL +rio_arg = xr.open_rasterio("soy_yield_1975_2016_005_1prc.tif") +rio_arg = rio_arg.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +rio_arg['time'] = list(map(int, rio_arg.attrs['descriptions'])) +rio_arg.name = 'Yield' + +ds_rio_arg = rio_arg.to_dataset() +ds_rio_arg = ds_rio_arg.sel(lat = slice(0,-35), lon = slice(-73.98,-34.72)) # SLICE IT TO BECOME FASTER TO CHANGE RESOLUTION LATER ON CDO +# ds_rio = ds_rio.where(ds_rio['harvest_area'] > -1000) + +ds_rio_arg['lat'].attrs = {'units':'degrees_north'} +ds_rio_arg['lon'].attrs = {'units':'degrees_east'} +ds_rio_arg['time'].attrs = {'units':f'years since {rio_arg.time[0]-1}-01-01'} +ds_rio_arg = ds_rio_arg.sortby('lat') +ds_rio_arg = ds_rio_arg.sortby('lon') +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio_arg.to_netcdf("soy_yield_1975_2016_005_1prc.nc") + +plot_2d_am_map(ds_rio_arg['Yield'].mean('time')) +plot_2d_am_map(ds_rio_arg['Yield'].sel(time = 1979)) + +#%% TRANSFORM usda + +# Load Geotiff +os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/Paper_drought/data') + +rio = xr.open_rasterio("soy_yields_US_all_1980_2020_05x05.tif") +rio = rio.rename({'y': 'lat', 'x': 'lon', 'band':'time'}) +rio['time'] = list(map(int, rio.attrs['descriptions'])) +rio.name = 'usda_yield' + +ds_rio = rio.to_dataset() +# ds_rio = ds_rio.where(ds_rio['harvest_area'] > -1000) +ds_rio['time'].attrs = {'units':'years since 1980-01-01'} +ds_rio['lat'].attrs = {'units':'degrees_north'} +ds_rio['lon'].attrs = {'units':'degrees_east'} +# ds_rio = ds_rio.mean('time') +# ds_rio_2 = ds_rio.sel(time = slice('1980', '2016')) +ds_rio.to_netcdf("soy_yields_US_all_1980_2020_05x05_02.nc") + +plot_2d_us_map(ds_rio['usda_yield'].mean('time')) + +df_rio = rio.to_dataframe().dropna() +df_rio_mean = df_rio.groupby('time').mean(...) + +#%% + +# TEST + +# cal_test = pd.read_csv('calendar_soybeans/mirca/CELL_SPECIFIC_CROPPING_CALENDARS_30MN.txt', sep = "\t") +# cal_test_sub = cal_test[(cal_test['crop'] == 34) & (cal_test['subcrop'] == 1)] +# cal_test_sub = cal_test_sub[['lat', 'lon', 'start', 'end']] +# cal_test_sub.set_index(['lat', 'lon'], inplace=True) + +# DS_cal_test_sub = cal_test_sub.to_xarray() +# DS_cal_test_sub.to_netcdf('mirca2000_soy_calendar.nc') + +# # Check for duplicates +# cal_test_sub[cal_test_sub.index.duplicated()] +# cal_test.query('lat==-39.25 & lon==-63.25 & crop==34') + +# plot_2d_us_map(DS_cal_test_sub['start']) +# plot_2d_us_map(DS_cal_test_sub['end']) + +#%% Load Geotiff +os.chdir('C:/Users/morenodu/OneDrive - Stichting Deltares/Documents/PhD/paper_hybrid_agri/data') +# TEST +DS_y_planting_flach = xr.open_dataset("soy_planting_area_1981_2019.nc", decode_times=True).sel(time = slice('1980', '2016')) +DS_y_harvested_flach = xr.open_dataset("soy_harvest_area_1981_2019.nc", decode_times=True).sel(time = slice('1980', '2016')) +# plot_2d_map(DS_y_planting_flach.Planted_area.mean('time')) +# plot_2d_map(DS_y_harvested_flach.Harvested_area.mean('time')) + +test = DS_y_planting_flach.Planted_area.sel(lon=slice(-55.12,-55.1), lat=slice(-13.63,-13.65), time=slice(1980,2014)) +test_b = test.sel(time=slice(2013,2014)) +test_c = test.sel(lon=slice(-55.11,-55.1), lat=slice(-13.64,-13.65)) +# Create a mock test where two cells are np.nan +counter_test = test_b.copy() +counter_test.values[0,0] = np.nan + +xr.where( (test_b > 0 ) & (np.isnan(counter_test) == True ), 1, 0 ) + +test_2=DS_y_planting_flach.Planted_area.sel(lon=slice(-74.25, -74.23), lat=slice(6.246, 6.237), time=slice(2014,2014)) + + +def earth_radius(lat): + ''' + calculate radius of Earth assuming oblate spheroid defined by WGS84 + + Input + --------- + lat: vector or latitudes in degrees + + Output + ---------- + r: vector of radius in meters + + Notes + ----------- + WGS84: https://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Chapter%203.pdf + ''' + from numpy import deg2rad, sin, cos + + # define oblate spheroid from WGS84 + a = 6378137 + b = 6356752.3142 + e2 = 1 - (b**2/a**2) + + # convert from geodecic to geocentric + # see equation 3-110 in WGS84 + lat = deg2rad(lat) + lat_gc = np.arctan( (1-e2)*np.tan(lat) ) + + # radius equation + # see equation 3-107 in WGS84 + r = ( (a * (1 - e2)**0.5) / (1 - (e2 * np.cos(lat_gc)**2))**0.5 ) + + return r + + +def area_grid(lat, lon): + """ + Calculate the area of each grid cell + Area is in square meters + + Input + ----------- + lat: vector of latitude in degrees + lon: vector of longitude in degrees + + Output + ----------- + area: grid-cell area in square-meters with dimensions, [lat,lon] + + Notes + ----------- + Based on the function in + https://github.com/chadagreene/CDT/blob/master/cdt/cdtarea.m + """ + from numpy import meshgrid, deg2rad, gradient, cos + from xarray import DataArray + + xlon, ylat = meshgrid(lon, lat) + R = earth_radius(ylat) + + dlat = deg2rad(gradient(ylat, axis=0)) + dlon = deg2rad(gradient(xlon, axis=1)) + + dy = abs(dlat) * R + dx = dlon * R * cos(deg2rad(ylat)) + + area = dy * dx + + xda = DataArray( + area, + dims=["lat", "lon"], + coords={"lat": lat, "lon": lon}, + attrs={ + "long_name": "area_per_pixel", + "description": "area per pixel", + "units": "m^2", + }, + ) + return xda + + +def fraction_mask(dataArray, percentage_limit): + + da_area = area_grid(dataArray['lat'], dataArray['lon']) + time_da = dataArray.time + + # Expand the mask to the entire time period + da_area_time = da_area.expand_dims(time=time_da) + + #Check if all values along time are equal -> True expected + print(np.all(np.diff(da_area_time, axis=0) == 0)) + + # Final mask with the fraction of harvested area by total grid area + area_frac = dataArray/da_area_time + + # Final mask selecting only grid cells above 5% of harvested area + da_clip = dataArray.where(area_frac > percentage_limit) + return da_clip + + +da_area = area_grid(DS_y_harvested_flach.Harvested_area['lat'], DS_y_harvested_flach.Harvested_area['lon']) +time_da = DS_y_harvested_flach.Harvested_area.time + +# Expand the mask to the entire time period +da_area_time = da_area.expand_dims(time=time_da) + +#Check if all values along time are equal -> True expected +print(np.all(np.diff(da_area_time, axis=0) == 0)) + +# Final mask with the fraction of harvested area by total grid area +area_frac = DS_y_harvested_flach.Harvested_area/da_area_time + +area_frac_perc = area_frac * 100 +area_frac_perc_nozero = area_frac_perc.where(area_frac_perc > 0) +area_frac_perc_nozero.plot.hist( bins = 100) +# sns.histplot(data=area_frac_perc) + +area_frac_perc_nozero_count_0 = xr.where(area_frac_perc_nozero > 0, 1, 0) +area_frac_perc_nozero_count_1 = xr.where(area_frac_perc_nozero > 1, 1, 0) + +#%% Different approaches for masks + +DS_y_harvested_flach_clip = fraction_mask(DS_y_harvested_flach.Harvested_area, 0.01) +DS_y_planting_flach_clip = fraction_mask(DS_y_planting_flach.Planted_area, 0.01) + +test_3 = DS_y_obs['Yield'].where((DS_y_harvested_flach_clip > 0) ) + +test_3.to_netcdf('soy_ibge_yield_1981_2019_1percent.nc') + +plot_2d_map(DS_y_harvested_flach_clip.mean('time')) +plot_2d_map(DS_y_planting_flach_clip.mean('time')) + + +DS_plant_count = xr.where(DS_y_planting_flach['Planted_area'] > 0 , 1, 0) +DS_harvest_count = xr.where(DS_y_harvested_flach['Harvested_area'] > 0 , 1, 0) +DS_yield_count = xr.where(DS_y_obs['Yield'] > 0 , 1, 0) + +# If planted is > 0 and harvest is 0, then area = 0 +plant_harvest_zeros = xr.where( (DS_plant_count == 1) & (DS_harvest_count==0) , 1, 0 ) + +# IF planted is 0 and harvest is > 0, raise error +plant_harvest_error = xr.where( (DS_plant_count == 0) & (DS_harvest_count==1) , 1, 0 ) + + +DS_yield_plant_clip = DS_y_obs['Yield'].where( DS_plant_count == 1 ) +DS_yield_harvest_clip = DS_y_obs['Yield'].where( DS_harvest_count == 1) + +DS_yield_error_harv = xr.where( (DS_y_obs['Yield'] > 0) & (DS_harvest_count==0), 1, 0) # This should be zero +DS_yield_error_plant = xr.where( (DS_y_obs['Yield'] > 0) & (DS_plant_count==0), 1, 0) # This should be zero too + +DS_combo = xr.where( (DS_yield_plant_clip > 0) & (DS_yield_harvest_clip > 0), 1, 0) + + +DS_yield_clip_plant_harv = DS_y_obs['Yield'].where( DS_combo == 1 ) + + +plot_2d_map(plant_harvest_zeros.mean('time')) +plot_2d_map(plant_harvest_error.mean('time')) +plot_2d_map(DS_yield_clip_plant_harv.mean('time')) + + +plot_2d_map(DS_y_obs['Yield'].mean('time')) +plot_2d_map(DS_y_obs['Yield'].sel(time=2016)) + +# Testing for observed data following 1% cliping +test_3 = DS_y_obs['Yield'].where((DS_y_harvested_flach_clip > 0) ) + +plot_2d_map(test_3.mean('time')) + +# Consider 0 cases where harvest is NA and planting area exists +test_4 = xr.where( plant_harvest_zeros == 1, 0, DS_y_obs['Yield']) + +plot_2d_map(test_4.mean('time')) + +# test_4.to_netcdf('soy_ibge_yield_1981_2019_clip.nc') + + +#TODO: CHECK THIS +# Select only grid cells without NAs along time axis: +test_5 = xr.where( np.isnan(DS_y_obs).any(dim = 'time') == False, DS_y_obs, np.nan ) +# test_5.drop('spatial_ref') +# test_5 = test_4.where(np.isnan(test_4).any() == False) +plot_2d_map(test_5['Yield'].mean('time')) + +# Select only grid cells with at least 10 years of data: +test_6 = xr.where( (np.isnan(DS_y_obs['Yield']) == False ).sum(dim='time') >= 5, DS_y_obs, np.nan ) +plot_2d_map(test_6['Yield'].mean('time')) +if len(test_6.coords) >3 : + test_6=test_6.drop('spatial_ref') +test_6.to_netcdf('soy_ibge_yield_1981_2019_5.nc') + +# Final filter method + +# 1) Select only grid cells with > 1% Harvest area fractions +test_filter_1 = xr.where( ds_rio_2['harvested_area'] >= 0, DS_y_obs['Yield'], np.nan ) + +# 2) Planting area > 0 & Harvest area is NA: Yield = 0 +test_filter_2 = xr.where( plant_harvest_zeros == 1, 0, DS_y_obs['Yield']) + +# 3) Remove grid cells with values repeating more than 3 times +# test_filter_3 = TODO + +# 4) Select only grid cells with minimum of 5 years per grid cell +test_filter_4 = xr.where( (np.isnan(DS_y_obs['Yield']) == False ).sum(dim='time') >= 5, DS_y_obs['Yield'], np.nan ) + +plot_2d_map(DS_y_obs['Yield'].mean('time')) +plot_2d_map(test_filter_2.mean('time')) +plot_2d_map(plant_harvest_zeros.mean('time')) +test_filter_2 = test_filter_2.drop('spatial_ref') +test_filter_2.to_netcdf('soy_yield_1981_2019_1prc_harvplan.nc') + +#%% +# hectares +DSharvest_area_spam = xr.open_dataset("spam2010v1r0_global_harvested-area_soyb_r.nc", decode_times=True) + +DS_yield_error_harv = xr.where( (DS_y_obs['Yield'] > 0) & (DS_harvest_count==0), 1, 0) # This should be zero +DS_yield_error_harv = xr.where( (DS_yield_count == 0) & (DS_harvest_count > 0), 1, 0) # This should be zero + + +da_area = area_grid(DSharvest_area_spam.Band1['lat'], DSharvest_area_spam.Band1['lon']) + +#Check if all values along time are equal -> True expected +print(np.all(np.diff(da_area, axis=0) == 0)) + +# Final mask with the fraction of harvested area by total grid area +area_frac = DSharvest_area_spam.Band1/da_area * 10000 + +area_frac_perc = area_frac * 100 +area_frac_perc_nozero = area_frac_perc.where(area_frac_perc > 0) +area_frac_perc_nozero.plot.hist( bins = 100) +sns.histplot(data=area_frac_perc) + +plot_2d_map(DSharvest_area_spam.Band1) + + + +# (Yield per grid * Harvested area) / sum( harvested _ area) +# This should be equal to mean(yield per grid) +total_area = DS_y_harvested_flach['Harvested_area'].sum(("lon", "lat")) +yield_weighted = (DS_y_obs['Yield'] * DS_y_harvested_flach['Harvested_area']) / total_area + +# Check if it matches +weights = DS_y_harvested_flach['Harvested_area'] / total_area +yield_weighted_2 = weights * DS_y_obs['Yield'] +yield_weighted_mean_2 = yield_weighted_2.sum(("lon", "lat")) # Should be identical to #1 + + +yield_weighted_mean = yield_weighted.sum(("lon", "lat")) +yield_mean = DS_y_obs['Yield'].mean(("lon", "lat")) + +yield_weighted_mean.plot(label = 'weighted') +yield_mean.plot(label = 'non-weighted') +plt.legend() +yield_weighted_mean.name = 'Yield' +if len(yield_weighted_mean.coords) >1 : + yield_weighted_mean=yield_weighted_mean.drop('spatial_ref') +df_yield_weighted_mean = yield_weighted_mean.to_dataframe() + +# removal with a 2nd order based on the CO2 levels +coeff = np.polyfit(df_co2.values.ravel(), df_yield_weighted_mean.values, 1) +trend = np.polyval(coeff, df_co2.values.ravel()) +df_yield_weighted_mean_det = pd.DataFrame( df_yield_weighted_mean['Yield'] - trend, index = df_yield_weighted_mean.index, columns = df_yield_weighted_mean.columns) + df_yield_weighted_mean.mean() + + +##### +test_stack = DS_y_harvested_flach['Harvested_area'][:1,:2,:2].copy() +test_stack.values = [[[10, 10],[10, 10]]] +##### + + +#%% GEnerate tiff or netcdf files from the vectors we have for each area +import rasterio +from rasterio.features import rasterize +from rasterio.transform import from_bounds +from shapely.geometry import box, mapping +import json + + +from geocube.api.core import make_geocube + +input_geopackage = gpd.read_file('soy_harvest_area_arg_1978_2019.gpkg') + + +out_grid = make_geocube( + vector_data=input_geopackage,measurements=["anio","superficie_cosechada_ha"],datetime_measurements ='anio', + group_by="anio",geom=json.dumps(mapping(box(-180, -90, 180, 90))), resolution=(-0.5, 0.5), +) + +out_grid["column_name"].rio.to_raster("my_rasterized_column.tif") + +out_grid = out_grid.rename({'y': 'lat', 'x': 'lon', 'anio':'time'}) +out_grid['time'].attrs = {'units':'years since 1978-01-01'} +out_grid['lat'].attrs = {'units':'degrees_north'} +out_grid['lon'].attrs = {'units':'degrees_east'} + + +plot_2d_am_map(out_grid['superficie_cosechada_ha'].sel(time = 2006)) + +out_grid['superficie_cosechada_ha'].sel(time = 2000).sum() / 10**6 diff --git a/code/wp3_sim.py b/code/wp3_sim.py index 9a6f9fd..250a17a 100644 --- a/code/wp3_sim.py +++ b/code/wp3_sim.py @@ -12,7 +12,7 @@ import pandas as pd import cartopy import cartopy.crs as ccrs -import matplotlib as mplf +import matplotlib as mpl mpl.rcParams['axes.spines.right'] = False mpl.rcParams['axes.spines.top'] = False