library(sf)
library(ggplot2)
library(rgdal)
library(dplyr)
library(gdalUtilities)
library(RColorBrewer)
library(ggnewscale)
library(paletteer)
library(glue)
#Referencias en
##https://jayrobwilliams.com/posts/2020/09/spatial-sql


# Elección de las zonas del área de México
# https://www.keene.edu/campus/maps/tool/
# Se guardan las coordenadas en archivo de texto y se convierten en un objeto espacial y con cierta proyección.
sf::sf_use_s2(FALSE)

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")

limitesbox2<-read.csv(paste0(rutacam,"2limites.txt"),header=TRUE, sep=",",encoding="latin")
names(limitesbox2)<-c("long","lat")

#Creación de zonas epaciales y su proyección

limitesbox$region <- "a"
limitesbox2$region <- "b"

dattot <- rbind(limitesbox,limitesbox2)

#Polígonos y proyección
sps <- dattot %>%
  st_as_sf(coords = c("long", "lat")) %>%
  group_by(region) %>%
  summarize(geometry = st_union(geometry)) %>%
  st_convex_hull()


st_crs(sps) <- CRS("+proj=longlat +datum=WGS84 +no_defs")

#Texto para consulta SQL
bb1 <- sps[1,] %>% 
  pull(geometry) %>% 
  st_as_text()

bb2 <- sps[2,] %>% 
  pull(geometry) %>% 
  st_as_text()

Conversión a archivos gpkg. Función :ogr2ogr library()

nomfile<-‘D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/mexico-latest-free.shp/gis_osm_roads_free_1.shp’ nomfil<-“gis_osm_roads_free_1” target_file = gsub(“.shp”, “.gpkg”, basename(nomfil)) ogr2ogr(src_datasource_name = nomfile, dst_datasource_name = target_file, f = “GPKG”)

Lectura del archivo con especificaciones de filtrado en consulta SQL

nomfil<-"gis_osm_roads_free_1"

a<-glue("SELECT * FROM ", nomfil, " WHERE st_intersects(geom, st_polygonfromtext(")
c<-paste0(a,"'",bb1,"' )) OR st_intersects(geom, st_polygonfromtext(","'",bb2,"' )) ")
recorteosmlayer<-read_sf(nomfil, query = c)

#plot(recorteosmlayer)
nomfil<-"gis_osm_waterways_free_1"
a<-glue("SELECT * FROM ", nomfil, " WHERE st_intersects(geom, st_polygonfromtext(")
#c<-paste0(a,"'",bb,"' ))")
c<-paste0(a,"'",bb1,"' )) OR st_intersects(geom, st_polygonfromtext(","'",bb2,"' )) ")
recorteosmlayer2<-read_sf(nomfil, query = c)

nomfil<-"gis_osm_railways_free_1"
a<-glue("SELECT * FROM ", nomfil, " WHERE st_intersects(geom, st_polygonfromtext(")
#c<-paste0(a,"'",bb,"' ))")
c<-paste0(a,"'",bb1,"' )) OR st_intersects(geom, st_polygonfromtext(","'",bb2,"' )) ")
recorteosmlayer3<-read_sf(nomfil, query = c)

nomfil<-"red_vial"
a<-glue("SELECT * FROM ", nomfil, " WHERE st_intersects(geom, st_polygonfromtext(")
#c<-paste0(a,"'",bb,"' ))")
c<-paste0(a,"'",bb1,"' )) OR st_intersects(geom, st_polygonfromtext(","'",bb2,"' )) ")
recorteosmlayer4<-read_sf(nomfil, query = c)

nomfil<-"gis_osm_pois_a_free_1"
a<-glue("SELECT * FROM ", nomfil, " WHERE st_intersects(geom, st_polygonfromtext(")
#c<-paste0(a,"'",bb,"' ))")
c<-paste0(a,"'",bb1,"' )) OR st_intersects(geom, st_polygonfromtext(","'",bb2,"' )) ")
recorteosmlayer5<-read_sf(nomfil, query = c)
#Gráfica de mapa
# librería ggnewscale para cambiar de paletas de colores

mapamx <- ggplot(data = recorteosmlayer,aes(color=fclass))          + geom_sf(show.legend = FALSE)+
  labs(title = "", x = "", y = "")+
   scale_colour_manual(values =via_paletacoltor)+
  ggnewscale::new_scale_color() +
  geom_sf(data=recorteosmlayer2,aes(color=fclass),show.legend = FALSE)+
  labs(title = "", x = "", y = "")+
   scale_colour_manual(values =via_paletagua)+
  ggnewscale::new_scale_color() +
  geom_sf(data=recorteosmlayer3,aes(color=fclass,linetype=factor(lintip),size=fclass),show.legend = FALSE)+
  labs(title = "", x = "", y = "")+
  scale_colour_manual(values =via_transporte)+
  #scale_colour_brewer(palette = "Greys")+
  scale_linetype_manual(values = LINES)+
  scale_size_manual(values=via_transporte2)+
  ggnewscale::new_scale_color() +
  geom_sf(data=recorteosmlayer4,aes(color=fclass),show.legend = FALSE)+
  labs(title = "", x = "", y = "")+
#  scale_colour_manual(values =via_transporte)+
 scale_colour_brewer(palette = "Greys")+
  ggnewscale::new_scale_color() +
  geom_sf(data=recorteosmlayer5,aes(color=fclass),show.legend = FALSE)+
  labs(title = "", x = "", y = "")+
#  scale_colour_manual(values =via_transporte)+
 scale_colour_brewer(palette = "Purples")+
 theme(
    panel.background = element_rect(fill = "transparent"), # bg of the panel
    plot.background = element_rect(fill = "white", color = NA), # bg of the plot
    panel.grid.major = element_blank(), # get rid of major grid
    panel.grid.minor = element_blank(), # get rid of minor grid
    legend.background = element_rect(fill = "transparent"), # get rid of legend bg
    legend.box.background = element_rect(fill = "transparent"), # get rid of legend panel bg
    legend.key = element_rect(fill = "transparent", colour = NA), # get rid of key legend fill, and of the surrounding
    axis.line = element_line(colour = "transparent"),
    axis.text.x=element_blank(), axis.ticks.x = element_blank(),
  axis.text.y = element_blank(), axis.ticks.y = element_blank())+ coord_sf()

mapamx