Qué necesito

Paquetes de R

 # carga los paquetes requeridos
x <- c("sf","ggplot2","rgdal","dplyr","raster","OpenStreetMap","tmap","tmaptools","broom","tidyverse","RStoolbox")

lapply(x, library, character.only = TRUE)

Mapa de la zona del contexto

Señalar los puntos de las líneas que forman el polígono

https://www.keene.edu/campus/maps/tool/

memory.limit(size=1000000)

rutacam<-"D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/"
limitesbox<-read.csv(paste0(rutacam,"limites.txt"),header=TRUE, sep=",",encoding="latin")
names(limitesbox)<-c("long","lat")
p = Polygon(limitesbox)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))
sp::proj4string(sps) <- CRS("+proj=longlat +datum=WGS84 +no_defs")
crs(sps)
plot(sps)

sps_df <- tidy(sps)
ggplot() +
  geom_path(data = sps_df,aes(x = long, y = lat, group = group))+ 
  coord_quickmap()

Datos abiertos de OpenStreetMap

https://download.geofabrik.de/ https://download.geofabrik.de/north-america/mexico.html

rutamc<-"D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/mexico-latest-free.shp"
rutacam<-"D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/mexico-latest-free.shp/"
mexosmlayer <- readOGR(rutamc,"gis_osm_waterways_free_1", use_iconv = TRUE, encoding = "latin1")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Documents\Claudia\Midropbox\Investigacion y escritos\mapasnacionales\mexico-latest-free.shp", layer: "gis_osm_waterways_free_1"
## with 82019 features
## It has 5 fields
sp::proj4string(mexosmlayer) <- CRS("+proj=longlat +datum=WGS84 +no_defs")
crs(mexosmlayer)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
plot(mexosmlayer)

Recorte de las capas osm en la zona de contexto

recorteosmlayer = mexosmlayer[sps, ]
plot(recorteosmlayer)

crs(recorteosmlayer)

#transformación en dataframe
recorteosmlayer$id <-rownames(recorteosmlayer@data)
recorteosmlayer_df <- tidy(recorteosmlayer, region = "id")
recorteosmlayer_df <- left_join(recorteosmlayer_df,recorteosmlayer@data,by = "id")

ggplot() +
  geom_path(data = recorteosmlayer_df,aes(x = long, y = lat, group = group,color=fclass))+ 
  coord_quickmap()

Mapas raster INEGI (elevaciones)

https://www.inegi.org.mx/app/geo2/elevacionesmex/

Mapas estatales

memory.limit(size=1000000)
##GDALinfo("D:/Documents/Claudia/Midropbox/Investigacion y ##escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/Coahui##la_r15m.tif")

##COAH_RASTER <- 
##  brick("D:/Documents/Claudia/Midropbox/Investigacion y ##escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/Coahuila_r15m.tif")
###str(COAH_RASTER)
###RGB = brick("path/RGB.tif")
###https://gis.stackexchange.com/questions/403427/rgb-raster-##plot-error-in-r

##sp::proj4string(COAH_RASTER) <- CRS("+proj=longlat +datum=WGS84 +no_defs")
##crs(COAH_RASTER)
#############################################
#generate a list of input rasters ("grids")
#pattern = "*.tif$" - filters for main raster files only and skips any associated files (e.g. world files)
grids <- list.files("D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/" , pattern = "*.tif$")


#######
# Reading raster files
rutaras<-"D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/"
rst <- lapply(list.files(rutaras , pattern = "*.tif$",full.names = T), FUN = brick)

  
rst_ext <- lapply(rst, FUN =  extent)

# Calculating full extent
xmin_rst <- c(); xmax_rst <- c(); ymin_rst <- c(); ymax_rst <- c();
for (i in 1:length(rst_ext)) {
    xmin_rst <- c(xmin_rst, rst_ext[[i]]@xmin)
    ymin_rst <- c(ymin_rst, rst_ext[[i]]@ymin)
    xmax_rst <- c(xmax_rst, rst_ext[[i]]@xmax)
    ymax_rst <- c(ymax_rst, rst_ext[[i]]@ymax)
}
full_extent <- extent(min(xmin_rst), max(xmax_rst),
                      min(ymin_rst), max(ymax_rst))

# Creating raster from fucll extent and first rasters' CRS and resolution
bounding_rst <- raster(full_extent, 
                   crs = crs(rst[[1]]),
                   res = res(rst[[1]]))

# Resampling rasters to match attributes of the bounding raster
rst_resampled <- lapply(X = rst, FUN = function(x) {
    target_rst <- crop(bounding_rst, x)
    resample(x, target_rst, method="bilinear")
})

# Creating mosaic
rst_mosaic <- do.call("merge", c(rst_resampled, fun = mean))
#rst_mosaic <- do.call("mosaic", c(rst_resampled, fun = mean))
######
sp::proj4string(rst_mosaic) <- CRS("+proj=longlat +datum=WGS84 +no_defs")


#Recorte del raster en la zona de coontext
recorterasterinegi <- crop(rst_mosaic, extent(sps), snap="out")      

fr <- rasterize(sps, recorterasterinegi)   
    lr <- mask(x=recorterasterinegi, mask=fr)


plot(lr)
plot(recorteosmlayer,add=T)

#Gráfica en ggplot
zonaproyectada <- projectRaster(recorterasterinegi,
                               crs =crs(recorteosmlayer))

recorterasterinegi_df <- as.data.frame(lr, xy = TRUE)

ggplot() +
    geom_raster(data = recorterasterinegi_df , aes(x = x, y = y, fill = layer)) +
#x, y = y)) +  
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
  coord_quickmap()

plot(recorteosmlayer)

e <- extent(recorterasterinegi)

ggplot() +
  geom_raster(data = recorterasterinegi_df , aes(x = x, y = y, fill = layer)) +
  #x, y = y, fill = Coahuila_r15m)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
  xlim(e[1:2]) + ylim(e[3:4]) +
  geom_path(data = recorteosmlayer_df,aes(x = long, y = lat, group = group,color=fclass))+
  coord_sf()

Mapas raster OSM

CBS_osm2 <- read_osm(sps, type="bing",zoom=14)
srtm_raster2 = as(CBS_osm2, "Raster") 
crs(srtm_raster2)
crs(recorteosmlayer)
zonaproyectada <- projectRaster(srtm_raster2,
                             crs = crs(recorteosmlayer))
plot(zonaproyectada)

#Recorte
e <- extent(zonaproyectada )
p <- ggplot() + ggRGB(img = zonaproyectada ,
                      r = 3,
                      g = 2,
                      b = 1,
                      stretch = 'hist',
                      ggLayer = T) + 
  xlim(e[1:2]) + ylim(e[3:4]) +
  geom_path(data = recorteosmlayer_df,aes(x = long, y = lat, group = group,color=fclass))+coord_sf()
p

Raster paquete Terra, maptiles

Obtener mosaicos de mapas desde distintos servidores

library(terra)
library(ggplot2)
library(tidyterra)
library(maptiles)
library(sf)
box <- st_as_sfc(sps)

tile <- get_tiles(box, provider = "OpenStreetMap", crop = TRUE, zoom = 16)
plot(tile)

#tile <- get_tiles(box, provider = "Wikimedia", crop = TRUE, zoom = 16)
#plot(tile)
#tile <- get_tiles(box, provider = "Esri.WorldStreetMap", crop = TRUE, zoom = 16)
#plot(tile)


#gráfica con ggplot
ggtile <- ggplot() +
  geom_spatraster_rgb(data = tile,maxcell = 500640)
ggtile

#spatvector = vect(recorteosmlayer)
#spatvector2 <- crop(spatvector, extent(sps))      


ggtile <- ggplot() +  geom_spatraster_rgb(data = tile,maxcell = 500640)+
  xlim(e[1:2]) + ylim(e[3:4]) +
  geom_path(data = recorteosmlayer_df,aes(x = long, y = lat, group = group,color=fclass))+coord_sf()
ggtile

#https://sesync-ci.github.io/blog/mapping-with-Mapbox.html