diff --git a/landform/landform.R b/landform/landform.R new file mode 100644 index 0000000..e1831aa --- /dev/null +++ b/landform/landform.R @@ -0,0 +1,105 @@ + + +# Load libraries ---------------------------------------------------------- +require(pacman) +pacman::p_load(RSAGA, terra, fs, ggrepel, sf, OpenStreetMap, colourpicker, rnaturalearth, rnaturalearthdata, tidyverse, rgeos, geodata, elevatr, gtools, glue, elevatr) + +g <- gc(reset = T) +rm(list = ls()) +options(scipen = 999, warn = -1) + +# Download ---------------------------------------------------------------- +wrld <- ne_countries(returnclass = 'sf', scale = 50, type = 'countries') +pan0 <- geodata::gadm(country = 'PAN', level = 0, path = './tmpr') +pan1 <- geodata::gadm(country = 'PAN', level = 1, path = './tmpr') +srtm <- geodata::elevation_30s(country = 'PAN', path = './tmpr') + +plot(srtm) + +# Extract by mask --------------------------------------------------------- +srtm <- terra::crop(srtm, pan0) +srtm <- terra::mask(srtm, pan0) + +dir_create('./tif') +terra::writeRaster(x = srtm, filename = './tif/srtm_raw.tif', overwrite = TRUE) + +srtm <- rast('./tif/srtm_raw.tif') +srtm <- terra::project(srtm, '+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +type=crs') + +terra::writeRaster(x = srtm, filename = './tif/srtm_prj.tif', overwrite = TRUE) + +srtm <- rast('./tif/srtm_prj.tif') + +# Rsaga ------------------------------------------------------------------- + +# Rsaga Environment +envr <- rsaga.env('C:/SAGA/saga-9.2.0_x64/saga-9.2.0_x64') + +# Rsaga check the parameters +rsaga.get.usage(lib = 'ta_morphometry', env = envr, module = 'TPI Based Landform Classification') + +# To execute TPI Based Landform Classification +rsl <- rsaga.geoprocessor( + lib = 'ta_morphometry', + module = 'TPI Based Landform Classification', + param = list(DEM = 'F:/yt/rsaga/landform/tif/srtm_prj.tif', + LANDFORMS = 'F:/yt/rsaga/landform/tif/landform.tif', + RADIUS_A_MIN = 0, + RADIUS_A_MAX = 100.000000, + RADIUS_B_MIN= 0.000000, + RADIUS_B_MAX= 1000.000000, + DW_WEIGHTING= 0), + env = envr) + +# To read the results +lndf <- terra::rast('./tif/landform.tif') +lndf <- terra::project(lndf, crs(pan0), method = 'near') + +# Landform labels +# Source: http://www.jennessent.com/downloads/tpi-poster-tnc_18x22.pdf +lbls <- tibble(value = 1:10, class = c('Streams', 'Midslope drainages', 'Upland drainages', 'Valleys', 'Plains', 'Open Slopes', 'Upper slopes', 'Local ridges', 'Midslope ridges', 'High ridges')) + +# Coordinates ------------------------------------------------------------ +crds <- pan1 %>% terra::centroids() %>% terra::crds() %>% as_tibble() %>% mutate(name_1 = pan1$NAME_1) + +# To make the map --------------------------------------------------------- +lndf.tble <- terra::as.data.frame(lndf, xy = T) %>% + as_tibble() %>% + setNames(c('x', 'y', 'value')) %>% + inner_join(., lbls, by = 'value') + +# Open Street Maps +ext(pan0) +LAT1 = 7.20235900000011 ; LAT2 = 9.6473610000001 +LON1 = -83.0518869819999 ; LON2 = -77.1292800909999 +map <- openmap(c(LAT2,LON1), c(LAT1,LON2), zoom = NULL, type = c("osm", "stamen-toner", "stamen-terrain","stamen-watercolor", "esri","esri-topo", 'esri-physical', 'esri-shaded')[8], mergeTiles = TRUE) +map.latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") + +gmap <- autoplot(map.latlon) + + geom_sf(data = wrld, fill = NA, col = 'grey30', inherit.aes = FALSE) + + geom_tile(data = lndf.tble, aes(x = x, y = y, fill = class)) + + scale_fill_manual(values = c('#BCBFA3', '#A5C2B2', '#BF793F', '#327D3C')) + + geom_sf(data = st_as_sf(pan0), fill = NA, col = 'grey40', inherit.aes = FALSE) + + geom_sf(data = st_as_sf(pan1), fill = NA, col = 'grey30', inherit.aes = FALSE) + + geom_text_repel(data = crds, aes(x = x, y = y, label = name_1), family = 'Gill Sans MT', size = 2, bg.color = 'white', bg.r = 0.25) + + ggtitle(label = 'Landform classification / Topograhic Position Index', + subtitle = 'Panama') + + coord_sf(xlim = ext(pan0)[1:2], ylim = ext(pan0)[3:4]) + + labs(x = 'Lon', y = 'Lat', fill = 'Landform', caption = 'CSI - SRTM') + + theme_minimal() + + theme(legend.position = 'bottom', + plot.title = element_text(face = 'bold', hjust = 0.5), + plot.subtitle = element_text(face = 'bold', hjust = 0.5), + text = element_text(family = 'Gill Sans MT'), + axis.text.y = element_text(angle = 90, hjust = 0.5)) + + guides(fill = guide_legend( + title.position = 'top', + title.hjust = 0.5, + nrow = 1 + )) + +gmap + +ggsave(plot = gmap, filename = './map_landform_pan.png', units = 'in', width = 12, height = 6, dpi = 300) + + \ No newline at end of file diff --git a/landform/landform.Rproj b/landform/landform.Rproj new file mode 100644 index 0000000..3af27f6 --- /dev/null +++ b/landform/landform.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/landform/map_landform_pan.png b/landform/map_landform_pan.png new file mode 100644 index 0000000..6f212a4 Binary files /dev/null and b/landform/map_landform_pan.png differ diff --git a/landform/tif/landform.tif b/landform/tif/landform.tif new file mode 100644 index 0000000..73c4cfe Binary files /dev/null and b/landform/tif/landform.tif differ diff --git a/landform/tif/srtm_prj.tif b/landform/tif/srtm_prj.tif new file mode 100644 index 0000000..b2c537b Binary files /dev/null and b/landform/tif/srtm_prj.tif differ diff --git a/landform/tif/srtm_raw.tif b/landform/tif/srtm_raw.tif new file mode 100644 index 0000000..a2fa2fe Binary files /dev/null and b/landform/tif/srtm_raw.tif differ diff --git a/landform/tmpr/PAN_elv_msk.tif b/landform/tmpr/PAN_elv_msk.tif new file mode 100644 index 0000000..5728715 Binary files /dev/null and b/landform/tmpr/PAN_elv_msk.tif differ diff --git a/landform/tmpr/gadm/gadm41_PAN_0_pk.rds b/landform/tmpr/gadm/gadm41_PAN_0_pk.rds new file mode 100644 index 0000000..12cdd88 Binary files /dev/null and b/landform/tmpr/gadm/gadm41_PAN_0_pk.rds differ diff --git a/landform/tmpr/gadm/gadm41_PAN_1_pk.rds b/landform/tmpr/gadm/gadm41_PAN_1_pk.rds new file mode 100644 index 0000000..84838ae Binary files /dev/null and b/landform/tmpr/gadm/gadm41_PAN_1_pk.rds differ diff --git a/osmMap_pro/toDownload.R b/osmMap_pro/toDownload.R new file mode 100644 index 0000000..7e4be9f --- /dev/null +++ b/osmMap_pro/toDownload.R @@ -0,0 +1,88 @@ + +## Fabio Alexander Castro Llanos +## Un geógrafo en YouTube + +# Load libraries ---------------------------------------------------------- +require(pacman) +pacman::p_load(sf, fs, glue, tidyverse, terra, osmdata, osrm) + +g <- gc(reset = T); rm(list = ls()); options(scipen = 999, warn = -1) + +# Function ---------------------------------------------------------------- +fetch_green <- function(key, values, boundary) { + + bbox <- st_bbox(boundary) + + q <- opq(bbox = bbox) |> + add_osm_feature(key = key, value = values) + + dat <- osmdata_sf(q) + + polys <- dat$osm_polygons + + if (!is.null(polys) && nrow(polys) > 0) { + polys <- st_transform(polys, st_crs(boundary)) + polys$key <- key + polys$value <- polys[[key]] + + polys_clipped <- st_intersection(polys, boundary) + return(polys_clipped) + } + + NULL + +} + +# To download ------------------------------------------------------------- +qbox <- opq(bbox = 'Cali') + +## Limits ---------------------------------------------- +lims <- qbox %>% add_osm_feature(key = "boundary", value = "administrative") %>% osmdata_sf() +lims <- lims$osm_multipolygons +cali <- lims %>% filter(name == 'Cali') + +### Levels Barrios + Comunas +lvls <- unique(lims$admin_level) +brrs <- lims %>% filter(admin_level == 9) +cmns <- lims %>% filter(admin_level == 8) + +## Road ------------------------------------------------- +road <- add_osm_feature(opq = qbox, key = 'highway') # See more about map feature: https://wiki.openstreetmap.org/wiki/Map_features +sort(available_features()) +road <- osmdata_sf(road) +road_line <- road$osm_lines +road_pnts <- road$osm_points +road_poly <- road$osm_polygons + +## Rivers ---------------------------------------------- +rvrs <- add_osm_feature(opq = qbox, key = 'waterway', value = 'river') +rvrs <- osmdata_sf(rvrs) +rvrs <- rvrs$osm_lines + +## Zonas verdes ---------------------------------------- +green_list_cali <- list(fetch_green(key = "leisure", values = c("park", "garden", "recreation_ground", "pitch", "nature_reserve"), boundary = cali), + fetch_green(key = "landuse", values = c("forest", "grass", "meadow"), boundary = cali), + fetch_green(key = "natural", values = c("wood", "heath", "scrub"), boundary = cali)) +saveRDS(object = green_list_cali, file = './rds/green_list_cali.rds') + +green_list_cali <- readRDS(file = './rds/green_list_cali.rds') +green_list_cali <- lapply(green_list_cali, function(x) {if(!is.null(x)) dplyr::select(x, osm_id, name, key, value, geometry)}) + +leis <- green_list_cali[[1]] +land <- green_list_cali[[2]] +natu <- green_list_cali[[3]] + +# To write the results --------------------------------------------------- +dir_create('./gpkg') +st_write(obj = cali, './gpkg/cali.gpkg') +st_write(obj = brrs, './gpkg/barrios.gpkg') +st_write(obj = cmns, './gpkg/comunas.gpkg') +st_write(obj = road_line, './gpkg/road_line.gpkg') +st_write(obj = road_poly, './gpkg/road_poly.gpkg') +st_write(obj = rvrs, './gpkg/rivers.gpkg') +st_write(obj = leis, './gpkg/leisure.gpkg') +st_write(obj = land, './gpkg/land.gpkg') +st_write(obj = natu, './gpkg/natu.gpkg') + +st_write(obj = lims, './gpkg/lims.gpkg') + diff --git a/trnd/climate graph.Rproj b/trnd/climate graph.Rproj new file mode 100644 index 0000000..3af27f6 --- /dev/null +++ b/trnd/climate graph.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/trnd/download climate.R b/trnd/download climate.R new file mode 100644 index 0000000..6a4591b --- /dev/null +++ b/trnd/download climate.R @@ -0,0 +1,131 @@ + +# Load libraries ---------------------------------------------------------- +require(pacman) +pacman::p_load(terra, sf, fs, tidyverse, glue) + +g <- gc(reset = T) +rm(list = ls()) +options(scipen = 999, warn = -1) + +# Load data --------------------------------------------------------------- +bsin <- vect('D:/data/spatial/igac/basins_col.gpkg') +dpto <- vect('D:/data/spatial/igac/mpios.gpkg') + +# Project the shapefile +bsin <- terra::project(bsin, crs(dpto)) + +# Intersection ------------------------------------------------------------ +bsin.dpto <- terra::intersect(bsin, dpto) +bsin.dpto +bsin.cali <- bsin.dpto[bsin.dpto$MPIO_CNMBR == 'CALI',] +bsin.cali$NOMSZH +bsin.cali$NOMSZH <- gsub('�', 'í', bsin.cali$NOMSZH) + +bsin.zone <- bsin.cali[bsin.cali$NOMSZH == 'Ríos Cali, río Lilí, río Melendez y río Canaveralejo',] +bsin.zone <- st_as_sf(bsin.zone) + +# To write +terra::writeVector(bsin.zone, './gpkg/rios_zone.gpkg', overwrite = T) +st_write(bsin.zone, './gpkg/bsin.gpkg') + +# To download ------------------------------------------------------------- + +# Parameters -------------------------------------------------------------- +root <- 'https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/NEX-GDDP-CMIP6' +ssps <- c('ssp245', 'ssp585') +vars <- c('pr', 'tasmax', 'tasmin') +gcms <- c('ACCESS-CM2', 'ACCESS-ESM1-5', 'BCC-CSM2-MR', 'CanESM5', 'CESM2-WACCM', 'CESM2', 'CMCC-CM2-SR5', 'CMCC-ESM2', 'CNRM-ESM2', 'CNRM-ESM2', 'CNRM-CM6-1', 'CNRM-ESM2-1', 'EC-Earth3-Veg-LR') + +# Historic dataset downlaod ---------------------------------------------- +down.hist <- function(var, gcm, ab1, ab2){ + + # Proof + var <- 'tasmax' + gcm <- gcms[1] + ab1 <- 'r1i1p1f1' + ab2 <- 'gn' + + # Start + cat('To process: ', var, ' ', gcm, '\n') + urlw <- as.character(glue('{root}/{gcm}/historical/{ab1}/{var}/{var}_day_{gcm}_historical_{ab1}_{ab2}_{1989}.nc')) + dirs <- as.character(glue('./tif/climate/cmip6_daily/historical/{gcm}/{var}/{basename(urlw)}')) + dout <- unique(dirname(dirs)) + ifelse(!file.exists(dout), dir_create(dout), print('Directorio existe')) + + # To download + map(.x = 1:length(urlw), .f = function(i){ + + cat('>>> To start the process:', i, '\t') + + # TO select the url for download + url <- urlw[i] + out <- dirs[i] + + # To download + download.file(url = url, destfile = out, mode = 'wb') + + rst <- rast(out) + + # To extract by mask + rst <- rotate(rst) + rst <- crop(rst, bsin.zone) + + # To write the final raster + terra::writeRaster(x = rst, filename = glue('{dout}/bsin_{basename(out)}'), overwrite = T) + file.remove(out); rm(rst); gc(reset = T) + cat('Done!\n') + + }) + +} + +# Future dataset download ------------------------------------------------- + +down.ftre <- function(var, gcm, ssp, ab1, ab2){ + + # Proof + # var <- 'tasmax' + # gcm <- gcms[1] + # ssp <- 'ssp370' + # ab1 <- 'r1i1p1f1' + # ab2 <- 'gn' + + # Start + cat('To process: ', var, ' ', gcm, '\n') + + 'https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/NEX-GDDP-CMIP6/ACCESS-CM2/ssp370/r1i1p1f1/tasmax/tasmax_day_ACCESS-CM2_ssp370_r1i1p1f1_gn_2015.nc' # R + 'https://nex-gddp-cmip6.s3-us-west-2.amazonaws.com/NEX-GDDP-CMIP6/ACCESS-CM2/ssp370/r1i1p1f1/tasmax/tasmax_day_ACCESS-CM2_ssp370_r1i1p1f1_gn_2015.nc' # WEB + + urlw <- as.character(glue('{root}/{gcm}/{ssp}/{ab1}/{var}/{var}_day_{gcm}_{ssp}_{ab1}_{ab2}_{2015:2016}.nc')) # 20monkemonke + + dirs <- as.character(glue('./tif/climate/cmip6_daily/future/{gcm}/{var}/{basename(urlw)}')) + dout <- unique(dirname(dirs)) + ifelse(!file.exists(dout), dir_create(dout), print('Directorio existe')) + + # To download + map(.x = 1:length(urlw), .f = function(i){ + + cat('>>> To start the process:', i, '\t') + + # To select the url for download + url <- urlw[i] + out <- dirs[i] + + # To download + download.file(url = url, destfile = out, mode = 'wb') + + rst <- rast(out) + + # To extract by mask + rst <- rotate(rst) + rst <- crop(rst, bsin.zone) + + # To write the final raster + terra::writeRaster(x = rst, filename = glue('{dout}/bsin_{basename(out)}'), overwrite = T) + file.remove(out); rm(rst); gc(reset = T) + cat('Done!\n') + + }) + +} + diff --git a/trnd/extract values.R b/trnd/extract values.R new file mode 100644 index 0000000..52cc977 --- /dev/null +++ b/trnd/extract values.R @@ -0,0 +1,89 @@ + +# Load libraries ---------------------------------------------------------- +require(pacman) +pacman::p_load(terra, sf, fs, tidyverse, glue, exactextractr) + +g <- gc(reset = T) +rm(list = ls()) +options(scipen = 999, warn = -1) + +# Load data --------------------------------------------------------------- + +# Vector data +bsin.zone <- vect('./gpkg/bsin.gpkg') + +# Raster data +path <- dir_ls('./tif/climate/cmip6_daily') + +# Baseline data ----------------------------------------------------------- +fles.bsln <- dir_ls(dir_ls(dir_ls(grep('historical', path, value = T)))) +fles.bsln <- grep('.nc$', fles.bsln, value = T) +fles.bsln <- as.character(fles.bsln) + +# Future data ------------------------------------------------------------- +fles.ftre <- dir_ls(dir_ls(dir_ls(grep('future', path, value = T))), regexp = '.nc$') +fles.ftre <- grep('.nc$', fles.ftre, value = T) +fles.ftre <- as.character(fles.ftre) + +# Extract values ---------------------------------------------------------- + +# Function to use +extr.vles <- function(fles){ + + # fles <- fles.bsln + + cat('... Starting ...\n') + rslt <- map(.x = fles, .f = function(f){ + + cat('To process: ', basename(f), '\n') + nme <- basename(f) + yea <- str_sub(nme, start = nchar(nme) - 6, end = nchar(nme) -3) + rst <- rast(f) + rst <- terra::crop(rst, bsin.zone) + rst <- terra::mask(rst, bsin.zone) + vle <- terra::as.data.frame(rst, xy = F) + vle <- as_tibble(t(vle)) + vle <- mean(pull(vle)) - 275.15 + vle <- tibble(year = yea, value = vle) + return(vle) + + }) + + rslt <- bind_rows(rslt) + return(rslt) + +} + +# To apply the function +vles.bsln <- extr.vles(fles = fles.bsln) +vles.ftre <- extr.vles(fles = fles.ftre) + +vles.bsln <- mutate(vles.bsln, Periodo = 'Línea base') +vles.ftre <- mutate(vles.ftre, Periodo = 'Futuro') +vles <- rbind(vles.bsln, vles.ftre) +vles <- mutate(vles, year = as.numeric(year)) +vles <- mutate(vles, Periodo = factor(Periodo, levels = c('Línea base', 'Futuro'))) + +plot(vles.bsln$value, type = 'l') +plot(vles$value, type = 'l') + +# To make the graph ------------------------------------------------------- + + +gline <- ggplot(data = vles, aes(x = year, y = value, col = Periodo)) + + geom_line() + + scale_x_continuous(labels = seq(1980, 2100, 10), breaks = seq(1980, 2100, 10)) + + scale_y_continuous(labels = 23:30, breaks = 23:30) + + scale_color_manual(values = c('#0B3B0B', '#DAA801')) + + ggtitle(label = 'Comportamiento de la temperatura máxima histórica y futura', + subtitle = 'SSP 370 - GCM: ACCESS-CM2') + + labs(x = 'Año', y = 'Temperatura (°C)', col = '') + + theme_light() + + theme(plot.title = element_text(face = 'bold', size = 14, hjust = 0.5), + plot.subtitle = element_text(face = 'bold', size = 14, hjust = 0.5), + legend.key.width = unit(3, 'line'), + legend.key.height = unit(1.2, 'line'), + text = element_text(family = 'Segoe UI'), + legend.position = 'bottom') + +ggsave(plot = gline, filename = './png/graph_trend.png', units = 'in', width = 9, height = 7, dpi = 300) diff --git a/worldData/download data.R b/worldData/download data.R new file mode 100644 index 0000000..d89d286 --- /dev/null +++ b/worldData/download data.R @@ -0,0 +1,76 @@ + +# Load libraries ---------------------------------------------------------- +require(pacman) +p_load(ggplot2, tidyverse, readxl, WDI, xlsx, FAOSTAT) + +# World indicators -------------------------------------------------------- +clss <- WDIsearch() %>% as_tibble() +View(clss) + +# Vars to download -------------------------------------------------------- +lbls <- tibble(vars = c('agricultural land', 'food inflation', 'population growth', 'food imports', 'food exports', 'foreing invest'), + abbr = c('AC.LND.AGRI.K2', NA, 'SP.POP.GROW', 'TM.VAL.FOOD.ZS.UN', 'TX.VAL.FOOD.ZS.UN', 'BX.KLT.DINV.CD.WD'), + source = c('Worldbank', 'FAO', 'Worldbank', 'Worldbank', 'Worldbank', 'Worldbank')) + +# Agricultural land = % of agricultural land +# Food inflation = % +# Pop growth = (%) +# Food imports = Food imports (% of merchandise imports) +# Food exports = Food exports (% of merchandise exports) +# Foreign invest = Foreign direct investment, net inflows (BoP, current US$)/ + +# To download ------------------------------------------------------------- +abbs <- lbls$abbr +abbs <- abbs[!is.na(abbs)] + +tble <- WDI(country = "all", indicator = abbs, start = 2001, end = 2020, extra = FALSE, cache = NULL, latest = NULL, language = "en") +tble <- as_tibble(tble) +tble <- dplyr::select(tble, country, iso = iso3c , year, SP.POP.GROW:BX.KLT.DINV.CD.WD) +colnames(tble)[4:7] <- c('pop_growth', 'food_imports', 'food_exports', 'foreing_invest') + +# Foreing invest +frgn <- WDI(country = "all", indicator = 'BX.KLT.DINV.CD.WD', start = 2001, end = 2020, extra = FALSE, cache = NULL, latest = NULL, language = "en") +frgn <- as_tibble(frgn) +frgn <- dplyr::select(frgn, iso = iso3c, year, foreing = BX.KLT.DINV.CD.WD) + +# Agricultural land (% of land area) +agrc <- WDI(country = "all", indicator = 'AC.LND.AGRI.K2', start = 2001, end = 2020, extra = FALSE, cache = NULL, latest = NULL, language = "en") +# Source: https://data.worldbank.org/indicator/AG.LND.AGRI.ZS +agrc <- read_csv('https://raw.githubusercontent.com/fabiolexcastro/tutorials_youtube/master/worldData/downloaded/API_AG.LND.AGRI.ZS_DS2_en_csv_v2_4669757/API_AG.LND.AGRI.ZS_DS2_en_csv_v2_4669757.csv', skip = 4) +agrc <- gather(agrc, var, value, -1, -2, -3, -4) +unique(agrc$`Indicator Name`) +unique(agrc$`Indicator Code`) +colnames(agrc) <- c('country_name', 'country_code', 'indicator_name', 'indicator_code', 'year', 'value') +agrc <- filter(agrc, year %in% as.character(2001:2020)) +agrc <- dplyr::select(agrc, country_name, country_code, indicator_name, year, value) +agrc <- dplyr::select(agrc, ISO = country_code, year, value) + +agrc <- mutate(agrc, year = as.numeric(year)) +tble +tble <- inner_join(tble, agrc, by = c('iso' = 'ISO', 'year' = 'year')) + +# Food inflation # Source: https://www.fao.org/faostat/en/#data/CP +food <- read_csv('https://raw.githubusercontent.com/fabiolexcastro/tutorials_youtube/master/worldData/downloaded/FAOSTAT_data_en_11-25-2022.csv') +colnames(food)[3] <- 'ISO' +food <- dplyr::select(food, ISO, Year, Item, Months, Value) +food <- food %>% group_by(ISO, Year, Item) %>% dplyr::summarise(Value = mean(Value, na.rm = T)) %>% ungroup() +food <- dplyr::select(food, -Item) +food <- rename(food, food_inflation = Value) + +tble <- tble %>% dplyr::rename(agricultural_land = value) +tble <- full_join(tble, food, by = c('iso' = 'ISO', 'year' = 'Year')) +write.xlsx(tble, file = './tble/database_global.xlsx', row.names = FALSE) + +# Deforestation +source_deforestation <- 'https://www.globalforestwatch.org/dashboards/global/?category=summary&location=WyJnbG9iYWwiXQ%3D%3D&map=eyJjZW50ZXIiOnsibGF0IjoxLjkwODMzMjgwODg3ODExZS0xMiwibG5nIjoxMS45OTk5OTk5OTk5OTcyNDN9LCJ6b29tIjoyLjIyMjE5MTAyNTA4OTE0NywiZGF0YXNldHMiOlt7ImRhdGFzZXQiOiJwb2xpdGljYWwtYm91bmRhcmllcyIsImxheWVycyI6WyJkaXNwdXRlZC1wb2xpdGljYWwtYm91bmRhcmllcyIsInBvbGl0aWNhbC1ib3VuZGFyaWVzIl0sImJvdW5kYXJ5Ijp0cnVlLCJvcGFjaXR5IjoxLCJ2aXNpYmlsaXR5Ijp0cnVlfSx7ImRhdGFzZXQiOiJOZXQtQ2hhbmdlLVNUQUdJTkciLCJsYXllcnMiOlsiZm9yZXN0LW5ldC1jaGFuZ2UiXSwib3BhY2l0eSI6MSwidmlzaWJpbGl0eSI6dHJ1ZSwicGFyYW1zIjp7InZpc2liaWxpdHkiOnRydWUsImFkbV9sZXZlbCI6ImFkbTAifX1dfQ%3D%3D&showMap=true' +dfrs <- read_csv('https://raw.githubusercontent.com/fabiolexcastro/tutorials_youtube/master/worldData/downloaded/P%C3%A9rdida%20mundial%20de%20bosques%20primarios/treecover_loss__ha.csv') +dfrs <- dfrs[,1:3] +colnames(dfrs) <- c('iso', 'year', 'loss_ha') # loss_ha = cambio de cobertura boscosa + +# To join both tables ----------------------------------------------------- +dfrs +tble <- as_tibble(full_join(tble, dfrs, by = c('iso' = 'iso', 'year' = 'year'))) + +# To write the final table ------------------------------------------------ +write.xlsx(tble, './tble/dabase_global_total.xlsx') +wqwwq diff --git a/worldData/make maps.R b/worldData/make maps.R new file mode 100644 index 0000000..afa05ad --- /dev/null +++ b/worldData/make maps.R @@ -0,0 +1,123 @@ + +# Load libraries ---------------------------------------------------------- +require(pacman) +p_load(ggplot2, tidyverse, classInt, showtext, colourpicker, extrafont, cptcity, RColorBrewer, sf, readxl, plotly, WDI, xlsx, FAOSTAT, rnaturalearth, rnaturalearthdata) + +g <- gc(reset = TRUE); rm(list = ls()); options(scipen = 999) + +# Font -------------------------------------------------------------------- +font_add_google(family = 'Roboto Condensed', name = 'Roboto Condensed') +showtext_auto() + +# Load data --------------------------------------------------------------- +shpf <- ne_countries(returnclass = 'sf', scale = 50) +tble <- read_excel(path = './tble/dabase_global_total.xlsx')[,-1] + +# Loss ha map ------------------------------------------------------------ +# Deforestation + +dfrs <- tble %>% group_by(iso) %>% summarise(loss_ha = sum(loss_ha, na.rm = T)) %>% ungroup() +nrow(dfrs); nrow(shpf) +shpf_dfrs <- inner_join(shpf, dfrs, by = c('adm0_a3_us' = 'iso')) +shpf_dfrs <- dplyr::select(shpf_dfrs, adm0_a3_us, geounit, loss_ha, geometry) +shpf_dfrs <- shpf_dfrs %>% drop_na() + +# find_cpt(name = 'hult'); image(matrix(1:100), col = cpt("hult_gr45_hult")) + +# To classify the colors +clss <- classIntervals(var = pull(shpf_dfrs, loss_ha), n = 6, style = 'fisher') +clss <- clss$brks +clss <- round(clss, -5) +clss <- c(clss[1], 500000, 1000000, clss[2:length(clss)]) +hist(pull(shpf_dfrs, loss_ha)) + +lbls <- tibble(value = 1:8, inf = clss[1:8], sup = clss[2:9], interval = paste0(inf, '-', sup)) +shpf_dfrs <- mutate(shpf_dfrs, class_lossHa = findInterval(x = shpf_dfrs$loss_ha, vec = clss, all.inside = TRUE)) +shpf_dfrs <- full_join(shpf_dfrs, lbls, by = c('class_lossHa' = 'value')) +shpf_dfrs <- mutate(shpf_dfrs, interval = factor(interval, levels = lbls$interval)) + +g_dfrs <- ggplot() + + geom_sf(data = shpf_dfrs, aes(fill = interval), lwd = 0.25, col = 'grey50') + + # scale_fill_viridis_d() + + scale_fill_manual(values = brewer.pal(n = 8, name = 'YlOrBr')) + + # scale_fill_gradientn(colors = rev(cpt(n = 10, 'hult_gr51_hult'))) + + labs(x = '', y = '', fill = 'Cobertura de bosque pérdida (ha)', caption = 'Adaptado de U. Maryland') + + ggtitle(label = 'Área en hectáreas de bosque deforestado entre 2000 y 2020') + + theme_minimal() + + theme(legend.position = 'bottom', + panel.background = element_rect(fill = 'white'), # 2C619E + panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + # legend.key.width = unit(3, 'line'), + legend.text = element_text(family = 'Roboto Condensed', size = 20), + axis.text.x = element_blank(), + plot.title = element_text(family = 'Roboto Condensed', size = 36, hjust = 0.5, face = 'bold'), + axis.text.y = element_blank(), + plot.caption = element_text(family = 'Roboto Condensed', size = 20), + axis.title = element_text(family = 'Roboto Condensed', size = 24, face = 'bold'), + legend.title = element_text(family = 'Roboto Condensed', size = 30)) + +ggsave(plot = g_dfrs, filename = 'png/global_deforested.png', units = 'in', width = 9, height = 6, dpi = 300) + +# Population growth ------------------------------------------------------- +popu <- tble %>% group_by(country, iso) %>% dplyr::summarise(pop_growth = mean(pop_growth, na.rm = T)) %>% ungroup() %>% drop_na() +popu <- inner_join(shpf, popu, by = c('adm0_a3_us' = 'iso')) + +g_popu <- ggplot() + + geom_sf(data = popu, aes(fill = pop_growth), lwd = 0.25, col = 'grey50') + + # scale_fill_manual(values = brewer.pal(n = 8, name = 'YlOrBr')) + + scale_fill_gradientn(colors = brewer.pal(n = 9, name = 'YlGnBu')) + + labs(x = '', y = '', fill = 'Crecimiento poblacional (%)', caption = 'Adaptado del Banco Mundial') + + ggtitle(label = 'Crecimiento promedio de la población') + + theme_minimal() + + theme(legend.position = 'bottom', + panel.background = element_rect(fill = 'white'), # 2C619E + panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + legend.key.width = unit(3, 'line'), + legend.text = element_text(family = 'Roboto Condensed', size = 20), + axis.text.x = element_blank(), + plot.title = element_text(family = 'Roboto Condensed', size = 36, hjust = 0.5, face = 'bold'), + axis.text.y = element_blank(), + plot.caption = element_text(family = 'Roboto Condensed', size = 20), + axis.title = element_text(family = 'Roboto Condensed', size = 24, face = 'bold'), + legend.title = element_text(family = 'Roboto Condensed', size = 30)) + +ggsave(plot = g_popu, filename = 'png/population_growth.png', units = 'in', width = 9, height = 6, dpi = 300) + +# Food inflation ---------------------------------------------------------- +finf <- tble %>% group_by(country, iso) %>% dplyr::summarise(food_inflation = mean(food_inflation, na.rm = T)) %>% ungroup() +finf <- drop_na(finf) +finf <- inner_join(shpf, finf, by = c('adm0_a3_us' = 'iso')) + +# To classify the colors +clss <- classIntervals(var = pull(finf, food_inflation), n = 6, style = 'fisher') +clss <- clss$brks +clss <- round(clss, 0) +lbls <- tibble(value = 1:6, inf = clss[1:6], sup = clss[2:7], interval = paste0(inf, '-', sup)) +finf <- mutate(finf, class_foodInflation = findInterval(x = finf$food_inflation, vec = clss, all.inside = TRUE)) +finf <- full_join(finf, lbls, by = c('class_foodInflation' = 'value')) +finf <- mutate(finf, interval = factor(interval, levels = lbls$interval)) + +g_finf <- ggplot() + + geom_sf(data = finf, aes(fill = interval), lwd = 0.25, col = 'grey50') + + # scale_fill_manual(values = brewer.pal(n = 8, name = 'YlOrBr')) + + scale_fill_manual(values = brewer.pal(n = 7, name = 'YlOrRd')) + + labs(x = '', y = '', fill = 'Inflación (%)', caption = 'Adaptado de la FAO') + + ggtitle(label = 'Inflación en los alimentos - Promedio entre 2000 y 2020') + + theme_minimal() + + theme(legend.position = 'bottom', + panel.background = element_rect(fill = 'white'), # 2C619E + panel.grid.major = element_blank(), + panel.grid.minor = element_blank(), + legend.text = element_text(family = 'Roboto Condensed', size = 20), + axis.text.x = element_blank(), + plot.title = element_text(family = 'Roboto Condensed', size = 36, hjust = 0.5, face = 'bold'), + axis.text.y = element_blank(), + plot.caption = element_text(family = 'Roboto Condensed', size = 20), + axis.title = element_text(family = 'Roboto Condensed', size = 24, face = 'bold'), + legend.title = element_text(family = 'Roboto Condensed', size = 30)) + +ggsave(plot = g_finf, filename = './png/food_inflation.png', units = 'in', width = 9, height = 6, dpi = 300) + +# Task: make the other maps.... and share with me in fabioalexandercastro@gmail.com diff --git a/xgb/out/run_1/glb_run1.RDS b/xgb/out/run_1/glb_run1.RDS new file mode 100644 index 0000000..c9a8c2a Binary files /dev/null and b/xgb/out/run_1/glb_run1.RDS differ diff --git a/xgb/out/run_1/glb_run2.RDS b/xgb/out/run_1/glb_run2.RDS new file mode 100644 index 0000000..7032505 Binary files /dev/null and b/xgb/out/run_1/glb_run2.RDS differ diff --git a/xgb/out/run_1/glb_run3.RDS b/xgb/out/run_1/glb_run3.RDS new file mode 100644 index 0000000..621d8c2 Binary files /dev/null and b/xgb/out/run_1/glb_run3.RDS differ diff --git a/xgb/out/run_1/glb_run4.RDS b/xgb/out/run_1/glb_run4.RDS new file mode 100644 index 0000000..17e2765 Binary files /dev/null and b/xgb/out/run_1/glb_run4.RDS differ diff --git a/xgb/out/run_1/glb_run5.RDS b/xgb/out/run_1/glb_run5.RDS new file mode 100644 index 0000000..183f081 Binary files /dev/null and b/xgb/out/run_1/glb_run5.RDS differ diff --git a/xgb/out/run_1/run_time.csv b/xgb/out/run_1/run_time.csv new file mode 100644 index 0000000..b8171f8 --- /dev/null +++ b/xgb/out/run_1/run_time.csv @@ -0,0 +1,12 @@ +"x" +2022-12-06 19:18:01 +2022-12-06 19:18:02,2022-12-06 19:18:36,33.6338307857513,"glb",1 +2022-12-06 19:18:36,2022-12-06 19:19:02,26.333477973938,"glb",2 +2022-12-06 19:19:02,2022-12-06 19:19:43,40.262393951416,"glb",3 +2022-12-06 19:19:43,2022-12-06 19:20:44,1.01509658495585,"glb",4 +2022-12-06 19:20:44,2022-12-06 19:21:32,47.586893081665,"glb",5 +2022-12-06 19:22:13,2022-12-06 19:48:23,26.1646463314692,"glb",1 +2022-12-06 19:48:23,2022-12-06 20:06:37,18.2345642964045,"glb",2 +2022-12-06 20:06:38,2022-12-06 20:19:39,13.0287772814433,"glb",3 +2022-12-06 20:19:40,2022-12-06 20:34:19,14.6529148022334,"glb",4 +2022-12-06 20:34:20,2022-12-06 20:50:08,15.8075969179471,"glb",5 diff --git a/xgb/out/run_2/glb_run1.RDS b/xgb/out/run_2/glb_run1.RDS new file mode 100644 index 0000000..ed56994 Binary files /dev/null and b/xgb/out/run_2/glb_run1.RDS differ diff --git a/xgb/out/run_2/glb_run2.RDS b/xgb/out/run_2/glb_run2.RDS new file mode 100644 index 0000000..77e3b2d Binary files /dev/null and b/xgb/out/run_2/glb_run2.RDS differ diff --git a/xgb/out/run_2/glb_run3.RDS b/xgb/out/run_2/glb_run3.RDS new file mode 100644 index 0000000..f3d2075 Binary files /dev/null and b/xgb/out/run_2/glb_run3.RDS differ diff --git a/xgb/out/run_2/glb_run4.RDS b/xgb/out/run_2/glb_run4.RDS new file mode 100644 index 0000000..5e2e5fc Binary files /dev/null and b/xgb/out/run_2/glb_run4.RDS differ diff --git a/xgb/out/run_2/glb_run5.RDS b/xgb/out/run_2/glb_run5.RDS new file mode 100644 index 0000000..990194a Binary files /dev/null and b/xgb/out/run_2/glb_run5.RDS differ diff --git a/xgb/out/run_2/run_time.csv b/xgb/out/run_2/run_time.csv new file mode 100644 index 0000000..184b0af --- /dev/null +++ b/xgb/out/run_2/run_time.csv @@ -0,0 +1,7 @@ +"x" +2022-12-10 08:54:55 +2022-12-10 08:55:02,2022-12-10 09:20:04,25.0200304508209,"glb",1 +2022-12-10 09:20:04,2022-12-10 09:38:21,18.2784571846326,"glb",2 +2022-12-10 09:38:21,2022-12-10 09:51:34,13.2156548817952,"glb",3 +2022-12-10 09:51:35,2022-12-10 10:06:28,14.8816575169563,"glb",4 +2022-12-10 10:06:28,2022-12-10 10:21:55,15.4379549185435,"glb",5 diff --git a/xgb/out/run_3/run_time.csv b/xgb/out/run_3/run_time.csv new file mode 100644 index 0000000..46414c7 --- /dev/null +++ b/xgb/out/run_3/run_time.csv @@ -0,0 +1,2 @@ +"x" +2022-12-10 11:48:19 diff --git a/xgb/xgb analysis.R b/xgb/xgb analysis.R new file mode 100644 index 0000000..f8eec00 --- /dev/null +++ b/xgb/xgb analysis.R @@ -0,0 +1,208 @@ + +# Load libraries ---------------------------------------------------------- +require(pacman) +pacman::p_load(tidyverse, xgboost, readxl, hrbrthemes, caret, naniar, glue, fs, sf, rnaturalearthdata, rnaturalearth, rmapshaper, fs, tidytext) + +g <- gc(reset = TRUE) +rm(list = ls()) +options(scipen = 999, warn = -1) + +# Functions --------------------------------------------------------------- + +# Function to find a "good" xgb model for each cluster +xgb.model.function <- function(level, run, cv_folds, cv_repeats, min_rsamp, tlength){ + + fseed <- run ^ 3 #seed defined using variable run + + # start time + time_start <- Sys.time() + + # print current level and run + cat(crayon::red(paste("starting\n[group: ", level, "]\n[run: run", run, "]\n", sep=""))) + + # subset data according to selected level + if(level %in% contn) { + dbfao_level = dbfao %>% + filter(continent == level) + print(paste("continent == ", level, sep="")) + } else if(level == "glb") { + dbfao_level = dbfao + print(paste("no filter, taking all contries", sep="")) + } else stop("problem with group selection in xgb.model.function") + + #print heads of data to be used + print(summary(dbfao_level)) + print(dbfao_level %>% arrange(iso3)) + + # select relevant columns + data_model <- dbfao_level %>% + select(df_ha, pop_growth, food_imports, food_exports, foreing_invest, + agricultural_land, food_inflation) + + # create training and testing data set + set.seed(fseed) # set seed for reproducibility + indata <- createDataPartition(y = data_model$df_ha, p = 0.8)[[1]] # index for testing and training data + training <- data_model[indata,] + testing <- data_model[-indata,] + + # trainControl object for repeated cross-validation with random search + adaptControl <- trainControl(method = "adaptive_cv", + number = cv_folds, + repeats = cv_repeats, + adaptive = list(min = min_rsamp, + alpha = 0.05, + method = "gls", + complete = FALSE), + search = "random", + verbose=TRUE, + allowParallel = TRUE) + + # train model + set.seed(fseed) # set seed for reproducibility + xgb_model <- train(df_ha ~ ., data = training, + method = "xgbTree", + trControl = adaptControl, + metric = "RMSE", + tuneLength = tlength, + na.action = na.pass) + + # model summary + model_summary <- xgb_model$results %>% + arrange(RMSE) %>% + head(1) %>% + transmute(rsquared = Rsquared, rmse = RMSE, eta, max_depth, gamma, colsample_bytree, + min_child_weight, subsample, nrounds) + + # prepare testing data + x_test <- select(testing, -df_ha) + y_test <- testing$df_ha + + # apply test data to model + predicted <- predict(xgb_model, x_test, na.action=na.pass) + + # calcuate residuals + test_residuals <- y_test - predicted + + # calculation of rsquare + tss <- sum((y_test - mean(y_test))^2) #total sum of squares + rss <- sum(test_residuals^2) #residual sum of squares + test_rsquared <- 1 - (rss/tss) #rsquare + + # calculate root mean square error + test_rmse <- sqrt(mean(test_residuals^2)) + # calculate mean absolut error + test_mae <- mean(abs(test_residuals)) + + # model evaluation + model_eval <- list(x_test=x_test, y_test=y_test, predicted=predicted, test_residuals=test_residuals, + test_rsquared=test_rsquared, test_rmse=test_rmse, test_mae=test_mae) + + # variable importance + model_var_importance <- + xgb.importance(feature_names = dimnames(data_model[,c(-1)])[[2]], model = xgb_model$finalModel) %>% + transmute(Feature, Gain) %>% + pivot_wider(names_from = Feature, values_from = Gain, names_prefix = "imp_") %>% + transmute(imp_pop_growth = ifelse(exists("imp_pop_growth"), imp_pop_growth, 0), + imp_food_imports = ifelse(exists("imp_food_imports"), imp_food_imports, 0), + imp_food_exports = ifelse(exists("imp_food_exports"), imp_food_exports, 0), + imp_foreing_invest = ifelse(exists("imp_foreing_invest"), imp_foreing_invest, 0), + imp_agricultural_land = ifelse(exists("imp_agricultural_land"), imp_agricultural_land, 0), + imp_food_inflation = ifelse(exists("imp_food_inflation"), imp_food_inflation, 0) + ) + + # save model settings + model_settings <- data.frame(level=level, seed=fseed, cv_folds=cv_folds, cv_repeats= cv_repeats, + min_rsamp=min_rsamp, tlength=tlength) + + # prepare output + output <- list(xgb_model, model_settings, model_summary, model_eval, model_var_importance, training, testing) + names(output) <- c(paste("xgb_model", level, sep = "_"), + paste("model_settings", level, sep = "_"), + paste("model_summary", level, sep = "_"), + paste("model_evaluation", level, sep = "_"), + paste("model_var_importance", level, sep = "_"), + paste("data_training", level, sep = "_"), + paste("data_testing", level, sep = "_")) + + # stop time + time_end <- Sys.time() + time_run <- time_end-time_start + + # print run time + cat(crayon::green(paste("[model run time: ", round(time_run, 2), "]\n", sep=""))) + + # write time file + run_time <- data.frame(time_start, time_end, time_run, level, run) + + write.table(run_time, file.path(path_output, "run_time.csv"), + sep = ",", col.names = FALSE, append=TRUE, row.names = FALSE) + + #return xgb output + assign(paste("output", level, fseed, sep = "_"), output) + +} + +# function to iterate over levels (clusters, geographical regions, continents) +xgb.run.function <- function(level, run, cv_folds, cv_repeats, min_rsamp, tlength){ + group <- ifelse(identical(level, contn), "contn", + ifelse(identical(level, glb), "glb", + stop("problem with group selection in xgb.run.function"))) + + run_x <- lapply(level, xgb.model.function, run, cv_folds, cv_repeats, min_rsamp, tlength) + names(run_x) <- level + saveRDS(run_x, file=file.path(path_output, paste(group, "_run", run, ".RDS", sep=""))) +} + +# Parameters -------------------------------------------------------------- +cv_folds = 10; cv_repeats = 10; min_rsamp = 5; tlength = 100 # Full +# cv_folds = 2; cv_repeats = 3; min_rsamp = 5; tlength = 6 # Testing +nrun <- 5 + +# Load data --------------------------------------------------------------- +dbfao <- read_excel('tbl/dabase_global_total.xlsx')[,-1] +path_output <- './out/run_3' +varbl <- c('df_ha', 'gdp_perc', 'food_exports', 'exports_gdp', 'food_infl', + 'pop_growth', 'gdp_growth', 'foreing_invest', 'tmp_change') + +wrld <- ne_countries(returnclass = 'sf', scale = 50) %>% dplyr::select(sov_a3, admin, region_un, geometry) +wrld_tble <- as_tibble(st_drop_geometry(wrld)) +dbfao <- inner_join(dbfao, wrld_tble, by = c('iso' = 'sov_a3')) +colnames(dbfao)[2] <- 'iso3' +colnames(dbfao)[12] <- 'continent' +colnames(dbfao)[10] <- 'df_ha' + +# Data preparation -------------------------------------------------------- +contn <- dbfao %>% + select(iso3, continent) %>% + unique() %>% + group_by(continent) %>% + dplyr::summarise(n = n()) + +# Define variable for whole world (no filtering) +glb <- "glb" +dbfao <- filter(dbfao, year != 2021) + +vis_miss(dbfao) +tail(dbfao) +n_miss(dbfao) +miss_case_table(dbfao) +prop_miss_var(dbfao) +miss_var_summary(dbfao) + +dbfao[is.na(dbfao$df_ha),] %>% View() +dbfao <- dbfao[!is.na(dbfao$df_ha),] + +# Run model --------------------------------------------------------------- +write.table(Sys.time(), file.path(path_output, "run_time.csv"), sep = ",", row.names = FALSE) + +# run with level glb +nrun <- 5 +level <- glb +lapply(c(1:nrun), xgb.run.function, level = level, cv_folds = cv_folds, cv_repeats = cv_repeats, min_rsamp = min_rsamp, tlength = tlength) # Normalizado + +# run with lemvel contn +level <- contn +lapply(c(1:nrun), xgb.run.function, level = level, cv_folds = cv_folds, cv_repeats = cv_repeats, min_rsamp = min_rsamp, tlength = tlength) + +# save run file and save ending time +write.table(Sys.time(), file.path(path_output, "run_time.csv"), sep = ",", col.names = FALSE, append=TRUE, row.names = FALSE) diff --git a/xgb/xgb results.R b/xgb/xgb results.R new file mode 100644 index 0000000..505150a --- /dev/null +++ b/xgb/xgb results.R @@ -0,0 +1,113 @@ + +# Load libraries ---------------------------------------------------------- +require(pacman) +pacman::p_load(tidyverse, readxl, xgboost, hrbrthemes, caret, glue, fs, sf, showtext, rnaturalearthdata, rnaturalearth, rmapshaper, fs, tidytext) + +g <- gc(reset = TRUE) +rm(list = ls()) +options(scipen = 999, warn = -1) + +# Fonts ------------------------------------------------------------------- +font_add(family = "Bahnschrift", regular = "Bahnschrift.ttf") +showtext_auto() + +# List the results -------------------------------------------------------- +dbse <- suppressMessages(read_excel('./tbl/dabase_global_total.xlsx')[,-1]) +path_output <- './out/run_2' + +# Global results +glb <- dir_ls(path_output, regexp = 'RDS') %>% grep('glb', ., value = TRUE) %>% as.character() %>% purrr::map(., readRDS) + +# Check the varnames +vrs <- glb[[1]]$glb$data_training_glb %>% colnames() + +# Variabels +lbels <- tibble(abb = c('df_ha', 'rsquared', 'pop_growth', 'food_imports', 'food_exports', 'foreing_invest', + 'agricultural_land', 'food_inflation'), + nme = c('Deforestation', 'R squared', 'Population growth', 'Food imports', 'Food exports', 'Foreing invest', + 'Agricultural land', 'Food inflation'), + nombre = c('Deforestación', 'R^2', 'Crecimiento poblacional', 'Importación alimentos', 'Exportación alimentos', 'Inversión extranjera', + 'Superficie agrícola', 'Inflación de los alimentos'), + units = c('Ha', '%', '%', '%', '% ', '(US$)', '%', '%') + ) + +lbels + +# Importance +prc <- list() +for(i in 1:5){ + prc[[i]] <- glb[[i]]$glb$model_var_importance_glb +} +prc <- bind_rows(prc) %>% mutate(run = LETTERS[1:5]) %>% gather(variable, value, -run) %>% mutate(run = factor(run, levels = LETTERS[1:5])) +prc %>% spread(run, value) + +# RSquared +rsq <- tibble(run = LETTERS[1:5], rsquared = unlist(map(.x = 1:5, .f = function(i)glb[[i]]$glb$model_summary_glb[1,1]))) +rsq <- mutate(rsq, rsquared = round(rsquared, 2)) +rsq <- transmute(rsq, value = rsquared, run, variable = 'rsquared') +rsq <- dplyr::select(rsq, run, variable, value) + +# Join importance table and rsquared table into only one table +all <- rbind(prc, rsq) +unique(all$variable) +unique(lbels$abb) +all <- mutate(all, + variable = factor(variable, + levels = unique(all$variable))) +all <- all %>% group_by(variable) %>% dplyr::summarise(min = min(value), avg = mean(value), max = max(value)) +all <- mutate(all, variable = gsub('imp_', '', variable)) +lbels +fix(lbels) +unique(all$variable) +lbels +all <- inner_join(all, lbels, by = c('variable' = 'abb')) + +write.csv(all, 'tbl/results_run2.csv', row.names = FALSE) + +# To make the graph ------------------------------------------------------ +all +# avg <- all %>% +# group_by(variable, nme, nombre, units) %>% +# dplyr::summarise(min = min(value), +# avg = mean(value), +# max = max(value)) %>% +# ungroup() +filter(avg, variable != 'rsquared') %>% pull(value) %>% sum() +avg <- all + +# To order ------------------------------------------------------------------ +lvls <- pull(arrange(avg, desc(avg)), variable) +avg <- mutate(avg, variable = factor(variable, levels = lvls)) +avg <- avg %>% arrange(variable) + +gbar <- ggplot(data = avg, aes(x = variable, y = avg, fill = variable, col = variable)) + + geom_col(position = 'dodge') + + geom_errorbar(aes(x = variable, ymin = min, ymax = max, group = variable), position = position_dodge(width = 0.9), width = 0.5, color = '#252525') + + scale_x_discrete(labels = avg$nombre) + + scale_fill_manual(values = c('#969696', "#d53e4f", "#f46d43", "#fdae61", "#e6f598", "#abdda4", "#66c2a5", "#3288bd", '#EDC949'), + labels = avg$nombre + ) + + scale_color_manual(values = c('#969696', "#d53e4f", "#f46d43", "#fdae61", "#e6f598", "#abdda4", "#66c2a5", "#3288bd", '#EDC949')) + + # geom_text(aes(y = avg, label = round(avg, 2)), vjust = 0, size = 10, fontface = 'bold', position = position_dodge(width = 0.9), col = 'grey40') + + ggtitle(label = 'Importancia de variables socio-económicas en función de la deforestación a nivel de país') + + theme_ipsum_ps() + + labs(x = '', y = 'Importance', fill = '', caption = 'Modelo XGBoost') + + theme(axis.text.x = element_text(angle = 45, family = 'Bahnschrift', size = 28, vjust = 0.5), + axis.text = element_text(family = 'Bahnschrift', size = 28), + axis.text.y = element_text(family = 'Bahnschrift', size = 28, face = 'bold'), + axis.title.x = element_text(family = 'Bahnschrift', face = 'bold', size = 28), + axis.title.y = element_text(family = 'Bahnschrift', face = 'bold', size = 28), + plot.title = element_text(family = 'Bahnschrift', face = 'bold', size = 40, hjust = 0.5), + panel.grid.major = element_blank(), + plot.caption = element_text(family = 'Bahnschrift', size = 28), + legend.text = element_text(family = 'Bahnschrift', size = 28), + legend.title = element_text(family = 'Bahnschrift', face = 'bold', size = 28), + strip.text = element_text(family = 'Bahnschrift', face = 'bold', hjust = 0.5, size = 28), + legend.background = element_rect(color = NA), + legend.position = 'bottom') + + guides(fill = guide_legend(nrow = 2), + color = 'none') + +gbar + +ggsave(plot = gbar, filename = 'png/graph.png', units = 'in', width = 9, height = 7, dpi = 300) diff --git a/xgb/xgbBoost.Rproj b/xgb/xgbBoost.Rproj new file mode 100644 index 0000000..3af27f6 --- /dev/null +++ b/xgb/xgbBoost.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX