fabiolexcastro · fabiolexcastro · Dec 5, 2022 · Dec 11, 2022 · Oct 20, 2023 · Nov 28, 2023
diff --git a/landform/landform.R 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)
+
+
diff --git a/landform/landform.Rproj 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
diff --git a/landform/tif/landform.tif b/landform/tif/landform.tif
diff --git a/landform/tif/srtm_prj.tif b/landform/tif/srtm_prj.tif
diff --git a/landform/tif/srtm_raw.tif b/landform/tif/srtm_raw.tif
diff --git a/landform/tmpr/PAN_elv_msk.tif b/landform/tmpr/PAN_elv_msk.tif
diff --git a/landform/tmpr/gadm/gadm41_PAN_0_pk.rds b/landform/tmpr/gadm/gadm41_PAN_0_pk.rds
diff --git a/landform/tmpr/gadm/gadm41_PAN_1_pk.rds b/landform/tmpr/gadm/gadm41_PAN_1_pk.rds
diff --git a/osmMap_pro/toDownload.R 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
@@ -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
@@ -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')
+
+  })
+
+}
+