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 https://www.keene.edu/campus/maps/tool/ https://www.inegi.org.mx/app/geo2/elevacionesmex/

GDALinfo("D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/Coahuila_r15m.tif")
## rows        38429 
## columns     29642 
## bands       1 
## lower left origin.x        -103.96 
## lower left origin.y        24.54268 
## res.x       0.0001388889 
## res.y       0.0001388889 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=longlat +ellps=GRS80 +no_defs 
## file        D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/Coahuila_r15m.tif 
## apparent band summary:
##   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1  Int16           TRUE       32767        128        128
## apparent band statistics:
##   Bmin Bmax    Bmean      Bsd
## 1  122 3706 1085.598 499.0108
## Metadata:
## AREA_OR_POINT=Area 
## DataType=Generic
memory.limit(size=1000000)
## [1] 1e+06
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)
## 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]]]]
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]]]]
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)
## 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]]]]
recorteosmlayer = mexosmlayer[sps, ]
plot(recorteosmlayer)

crs(recorteosmlayer)
## 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]]]]
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()

recorterasterinegi <- crop(COAH_RASTER, extent(sps), snap="out")      

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

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

zonaproyectada <- projectRaster(recorterasterinegi,
                               crs =crs(recorteosmlayer))
recorterasterinegi_df <- as.data.frame(lr, xy = TRUE)
#str(recorterasterinegi_df)
ggplot() +
    geom_raster(data = recorterasterinegi_df , aes(x = x, y = y, fill = Coahuila_r15m)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) + 
  coord_quickmap()

ggplot() +
  geom_raster(data = recorterasterinegi_df , aes(x = x, y = y, fill = Coahuila_r15m)) +
     scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
  geom_path(data = recorteosmlayer_df,aes(x = long, y = lat, group = group,color=fclass))+ 
  coord_equal()

CBS_osm2 <- read_osm(sps, type="bing",zoom=18)
srtm_raster2 = as(CBS_osm2, "Raster") 
crs(srtm_raster2)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1
## +units=m +nadgrids=@null +wktext +no_defs 
## WKT2 2019 representation:
## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["unnamed",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]
crs(recorteosmlayer)
## 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]]]]
zonaproyectada <- projectRaster(srtm_raster2,
                             crs = crs(recorteosmlayer))
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

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

tile <- get_tiles(box, crop = TRUE, zoom = 18)
ggtile <- ggplot() +
  geom_spatraster_rgb(data = tile,maxcell = 500640)
ggtile

spatvector = vect(recorteosmlayer)


ggtile <- ggplot() +
  geom_spatraster_rgb(data = tile,maxcell = 501215)+
  geom_spatvector(data = spatvector)+coord_sf()
ggtile