Vegetación desde dos fuentes

Ubicación espacial desde el estado y municipio.

Mapa visual Identificación de estado y municipio https://www.inegi.org.mx/contenidos/productos/prod_serv/contenidos/espanol/bvinegi/productos/geografia/imagen_cartografica/1_2_000_000/794551132258.pdf

Ejemplo:Jalisco 14 Guadalajara 039

Una vez ubicado el o los estados y municipios la identificación del polígon utilizando la información de los asentamientos humanos del INEGI

Recuperación de archivo INEGI

Colonias https://www.inegi.org.mx/programas/dcah/

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"
infile <- "00as.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"00as"

a<-"SELECT * FROM \"00as\" WHERE (CVE_ENT = '05' AND CVE_MUN = '035') OR (CVE_ENT='10' AND CVE_MUN ='007') OR (CVE_ENT='10' AND CVE_MUN='012')"
#a<-"SELECT * FROM \"00as\" WHERE (CVE_ENT = '14') AND (CVE_MUN='039')"

`%notin%` <- Negate(`%in%`)
sf_use_s2(FALSE)

datageot<- read_sf(nomarchi,options = "ENCODING=latin1",query=a)

datageot= st_transform(datageot, crs = 4326)
    datageot<- datageot|> filter(TIPO %notin% c("NINGUNO","CIUDAD"))


lim<-st_bbox(datageot)
lim$ymax<-25.65

xlim <- c(lim$xmin,lim$xmax)
ylim <- c(lim$ymin,lim$ymax)
colors <- paletteer::paletteer_d("khroma::soil")

g <- ggplot() +
  geom_sf(data = datageot, aes(fill = TIPO),alpha=.7) +
  scale_fill_manual(values = colors)+

  labs(
    title = "Asentamientos humanos ") +
  theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +
  theme(
    legend.position = "bottom",
    legend.key.width = unit(3, "line"),
    legend.key.height = unit(0.5, "line"),
    plot.background = element_rect(fill = NA, color = NA),
    plot.title = element_text(size = 20,  hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(size = 16,  hjust = 0.5),
    plot.caption = element_text(size = 10,  hjust = 0.5, margin = margin(20, 0, 0, 0)))+
  coord_sf(xlim = xlim, ylim = ylim)
g

Global canopy height map for the year 2020 derived from Sentinel-2 and GEDI

https://www.research-collection.ethz.ch/handle/20.500.11850/609802 Se descarga la tile https://libdrive.ethz.ch/index.php/s/cO8or7iOe5dT2Rt primero arhvivo tile_index.html para elegir tile file:///C:/Users/cguer/Downloads/tile_index.html

y se localiza economics la carpeta 3deg_cogs https://libdrive.ethz.ch/index.php/s/cO8or7iOe5dT2Rt?path=%2F%2F3deg_cogs

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/mapaschallenger/data/"
nomfile<-"ETH_GlobalCanopyHeight_10m_2020_N24W105_Map.TIF"
#nomfile<-"ETH_GlobalCanopyHeight_10m_2020_N18W105_Map.TIF"

nocom<-paste0(ruta,nomfile)
ras <- terra::rast(nocom)
#Torreón ETH_GlobalCanopyHeight_10m_2020_N24W105_Map
####### proyeccion de rasters.....importante
#""
#st_crs(ras)
#crs(ras, describe = TRUE)
wgs84 = "EPSG:4326"
ras_wgs84 = terra::project(ras, wgs84)

|———|———|———|———|=========================================

#plot(ras_wgs84)
library(swatches)
library(cptcity)

maparaster_a <- crop(ras,datageot)
#plot(maparaster_a)



maparaster_a_df <- as.data.frame(maparaster_a, xy = T)
maparaster_a_df <- na.omit(maparaster_a_df)
colnames(maparaster_a_df)[3] <- "height"

    maparaster_a_df<- maparaster_a_df|> filter(y < 25.65)


p <- ggplot(maparaster_a_df) +
  geom_raster(aes(x = x,y = y,fill = height)) +
  #scale_fill_gradientn(name = "altura (m)",colors = texture,breaks = round(breaks, 0)) +
  # scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 26, "neota_flor_burnt_tree"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 101, "neota_flor_grass"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 101, "pj_7_christmattree1"))+
  scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 254, "ncl_WhiteGreen"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 254, "gps_GPS_Nature_Random_Greens"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 26, "es_emerald_dragon_es_emerald_dragon_25"))+
    geom_sf(data = datageot, fill=NA,color="navy",linewidth = 0.5) +

  theme_minimal() +
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "bottom",
    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)
p

Localización de vialidades

q <- opq(bbox = c(lim$xmin, lim$ymin, lim$xmax, lim$ymax))%>% 
  add_osm_feature(key = "highway") %>% 
  osmdata_sf()
qosm_sf <- q$osm_lines 
st_crs(qosm_sf) <-4326


qosm_sf <- qosm_sf|> filter(highway %in% c("primary","secondary","tertiary","residential"))




shapename <- read_sf('C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/mexico-latest-free.shp/gis_osm_waterways_free_1.shp')
#filtros de tipos de corrientes de agua
#mexosmlayer2_sub <- shapename[shapename$fclass %in% c("river","stream"), ]
mexosmlayer2_sub<-shapename[datageot, ] 


p <- ggplot(maparaster_a_df) +
  geom_raster(aes(x = x,y = y,fill = height)) +
  #scale_fill_gradientn(name = "altura (m)",colors = texture,breaks = round(breaks, 0)) +
  # scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 26, "neota_flor_burnt_tree"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 101, "neota_flor_grass"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 101, "pj_7_christmattree1"))+
  scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 254, "ncl_WhiteGreen"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 254, "gps_GPS_Nature_Random_Greens"))+
  #scale_fill_gradientn(name = "altura (m)",colours = cpt(n = 26, "es_emerald_dragon_es_emerald_dragon_25"))+
   #geom_sf(data = recorteosmlayer4, fill = NA,alpha=0.1,linewidth=0.05,color="grey") +
    #geom_sf(data = datageot, fill=NA,color="navy",linewidth = 0.5) +
      geom_sf(data = qosm_sf, color="grey",fill = NA, linewidth = 0.1) +
  geom_sf(data = mexosmlayer2_sub, color="blue",fill = "blue",alpha=3) +

  theme_minimal() +
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
   # legend.position = "bottom",
    legend.position = "none",

    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)
p

Índice de Vegetación de Diferencias Normalizadas (NDVI) Landsat 8,9 del periodo correspondiente al año(s) 2023 del índice: MX_013007

https://www.inegi.org.mx/investigacion/ndvi/#mapas

Para localizar el mosaico correspondiente de observación

Mosaicos del cubo de datos geoespaciales de México CDGM

https://www.inegi.org.mx/app/biblioteca/ficha.html?upc=998880013761

https://www.inegi.org.mx/contenidos/productos/prod_serv/contenidos/espanol/bvinegi/productos/nueva_estruc/889463908272.pdf

Banda 1 Mí nimo Banda 2 Media Banda 3 Máximo Banda 4 Desv Banda 5 Mediana

La capa 1: el valor del índice mínimo registrado en el periodo de estudio. · La capa 2: el valor promedio del índice en periodo de estudio. · La capa 3: el valor máximo del índice en el periodo de estudio. · La capa 4: la desviación estándar del índice en el periodo de estudio. · La capa 5: el valor mediano del índice en el periodo de estudio

library(raster)
library(viridis)
ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"

nomfile2<-"ndvi_2023_MX_013007.tif"
#nomfile<-"ndvi_2023_MX_012010.tif"
#nomfile<-"ndvi_2023_MX_014007.tif"
nomfile1<-"ndvi_2023_MX_013006.tif"

nocom<-paste0(ruta,nomfile1)
ras <- terra::rast(nocom)

nocom<-paste0(ruta,nomfile2)
ras2 <- terra::rast(nocom)

rastot<-merge(ras,ras2)

|———|———|———|———|=========================================

#st_crs(rastot)
#crs(rastot, describe = TRUE)
wgs84 = "EPSG:4326"
ras_wgs84 = terra::project(rastot, wgs84)

|———|———|———|———|=========================================

ras_wgs84_a <- crop(ras_wgs84,datageot)


plot(ras_wgs84_a)

#https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
#rss <- "/ncdffilename.nc"##  having 60 layers (it's a raster with one single variable for 60 years), the variable name is "cdd"
#names(rss) <- "cdd" ## not necessary
#r1 <- terra::rast(rss,subds = "cdd") ## select the variable name
#r <- r1[[1]] ## give me the first layer
plot(ras_wgs84_a[[3]]) #MAXIMO

ras_wgs84_df <- as.data.frame(ras_wgs84_a, xy = TRUE, na.rm = TRUE)
    ras_wgs84_df<- ras_wgs84_df|> filter(y < 25.65)

rownames(ras_wgs84_df) <- c()
colnames(ras_wgs84_df)[3] <- "minim"
colnames(ras_wgs84_df)[4] <- "media"
colnames(ras_wgs84_df)[5] <- "maxim"
colnames(ras_wgs84_df)[6] <- "stdim"
colnames(ras_wgs84_df)[7] <- "medim"
ras_wgs84_df$maxim2<-ifelse(ras_wgs84_df$maxim<0,0,ras_wgs84_df$maxim)
  #gain(data=ras_wgs84_a) <- 0.0001
ggplot() +
  geom_raster(
    data = ras_wgs84_df,
    aes(x = x, y = y, fill = maxim2)
  ) +
   scale_fill_gradientn(name = "NDVI",colours = cpt(n = 254, "ncl_WhiteGreen"))+

  geom_sf(data = datageot, inherit.aes = FALSE, fill = NA,color="white") +

  labs(
    title = "NDVI (Normalized Difference Vegetation Index) in la Laguna",
    subtitle = "01-06-2020",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()+
  coord_sf(xlim,ylim)

p <- ggplot(ras_wgs84_df) +
  geom_raster(aes(x = x,y = y,fill = maxim)) +
   scale_fill_gradientn(name = "NDVI ",colours = cpt(n = 254, "ncl_WhiteGreen"))+
  geom_sf(data = qosm_sf, color="grey",fill = NA, linewidth = 0.1) +
  geom_sf(data = mexosmlayer2_sub, color="blue",fill = "blue",alpha=3) +

  theme_minimal() +
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "bottom",
   #legend.position = "none",

    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)
p

nocom<-paste0(ruta,nomfile1)
#NDVI_raster = raster(nocom)

NDVI_raster = brick(nocom)
NDVI_raster <- projectRaster(NDVI_raster, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
NDVI_raster <- raster::mask(NDVI_raster, as_Spatial(datageot))
gain(NDVI_raster) <- 0.0001
NDVI_df <- as.data.frame(NDVI_raster, xy = TRUE, na.rm = TRUE)
rownames(NDVI_df) <- c()

nocom<-paste0(ruta,nomfile2)
#NDVI_raster2 = raster(nocom)

NDVI_raster2 =brick(nocom)
NDVI_raster2 <- projectRaster(NDVI_raster2, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
NDVI_raster2 <- raster::mask(NDVI_raster2, as_Spatial(datageot))
gain(NDVI_raster2) <- 0.0001
NDVI_df2 <- as.data.frame(NDVI_raster2, xy = TRUE, na.rm = TRUE)
rownames(NDVI_df2) <- c()


#rastot<-sum(NDVI_raster,NDVI_raster2)
#NDVI_raster<-rastot


# Visualising using ggplot2
ggplot() +
  geom_raster(
    data = NDVI_df,
    aes(x = x, y = y, fill = ndvi_2023_MX_013006_3)
  ) +
  geom_raster(
    data = NDVI_df2,
    aes(x = x, y = y, fill = ndvi_2023_MX_013007_3)
  ) +
  geom_sf(data = qosm_sf, inherit.aes = FALSE, fill = NA,color="grey",linewidth = 0.1) +
    scale_fill_gradientn(name = "NDVI",colours = cpt(n = 254, "ncl_WhiteGreen"))+

 # scale_fill_viridis(name = "NDVI") +
  labs(
    title = "NDVI (Normalized Difference Vegetation Index) in La Laguna",
    subtitle = "01-06-2020",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()+
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "bottom",
   #legend.position = "none",
        legend.key.width = unit(4, 'cm'), #change legend key width

    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)

Otras fuentes de información https://rspatialdata.github.io/vegetation.html#MODIS_(Moderate_Resolution_Imaging_Spectroradiometer)

Información de la capa de INEGI

SIA Servicios con Información complementaria de tipo Área (áreas Verdes, Camellones, glorietas) SIP Servicios con Información complementaria de tipo Puntual (Palacios Municipales o Ayudantías, Parques o Jardines)

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/portraits/torreon/"

infile <- "malla_lagunanivel10.shp"
nomarchi<-paste0(ruta,infile)
#st_write(recorteosmlayer4, nomarchi)
#write_sf(recorteosmlayer4, nomarchi)

recorteosmlayer4<-read_sf(nomarchi)


sf_use_s2(FALSE)
ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"
infile <- "05sia.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"05sia"
a<-"SELECT * FROM \"05sia\" WHERE CVE_ENT = '05' AND CVE_MUN='035'"

datasia<- read_sf(nomarchi,options = "ENCODING=latin1",query=a)
datasia= st_transform(datasia, crs = 4326)
#unique(datasia$TIPO)
datasia<- datasia|> filter(TIPO %in% c("Jardín","Área Verde","Parque"))
recortedatasia <-datasia[datageot, ] 

#datageo_geom <- st_geometry(datageot)
#polygon = st_cast(datageo_geom, "POLYGON") 
#tor<-polygon[2]



ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/conjunto_de_datos/"
infile <- "10sia.shp"
nomarchi<-as.filename(paste0(ruta,infile))
nomfil<-"10sia"
a<-"SELECT * FROM \"10sia\" WHERE CVE_MUN = '012' OR CVE_MUN='007'"
#datasiaD<- read_sf(nomarchi,options = "ENCODING=latin1")

datasiaD<- read_sf(nomarchi,options = "ENCODING=latin1",query=a)
datasiaD= st_transform(datasiaD, crs = 4326)
#unique(datasia$TIPO)
datasiaD<- datasiaD|> filter(TIPO %in% c("Jardín","Área Verde","Parque"))
recortedatasiaD <-datasiaD[datageot, ] 



g <- ggplot() +
 
 geom_raster(
    data = NDVI_df,
    aes(x = x, y = y, fill = ndvi_2023_MX_013006_3)
  ) +
  geom_raster(
    data = NDVI_df2,
    aes(x = x, y = y, fill = ndvi_2023_MX_013007_3)
  ) +
  geom_sf(data = qosm_sf, inherit.aes = FALSE, fill = NA,color="navy",linewidth = 0.1) +
    geom_sf(data = recorteosmlayer4, inherit.aes = FALSE, fill = NA,color="grey",linewidth = 0.1) +

  
    scale_fill_gradientn(name = "NDVI",colours = cpt(n = 254, "ncl_WhiteGreen"))+
 geom_sf(data = recortedatasia, color = "#A45A52") +
    geom_sf(data = recortedatasiaD, color = "#A45A52") +
    theme_minimal(
    base_family = "body_font",
    base_line_size = unit(0.1, "mm")
  ) +

  labs(
    title = "NDVI SIA INEGI en La Laguna",
    subtitle = "01-06-2020",
    caption = "Parques jardines INEGI color #A45A52",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()+
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "bottom",
   #legend.position = "none",
        legend.key.width = unit(4, 'cm'), #change legend key width

    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)

g

#ggsave("veglocal.svg", width = 1300, height = 1300, units = "mm", dpi = "retina",limitsize = FALSE)

Información de parques del INEGI (mejor la capa sia no usar)

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/portraits/torreon/"

infile <- "INEGI_areasVerdesParques_.csv"
nomarchi<-paste0(ruta,infile)

parquesinegi<-read.csv(nomarchi)
parquesinegi$geometry <- st_as_sfc(parquesinegi$GEOM,crs=4326)



g <- ggplot() +
  

  geom_sf(data = parquesinegi$geometry, inherit.aes = FALSE, fill = "navy",color="navy",linewidth = 0.1) +
geom_sf(data = recortedatasia, fill= "#A45A52",color = "#A45A52") +
  theme_minimal()+
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "bottom",
   #legend.position = "none",
        legend.key.width = unit(4, 'cm'), #change legend key width

    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)

g

Información de parques del INEGI capa SIA Y OSM FEATURES Y parques y jardines INEGI

#•  boundary=national_park - large, government-decreed area for human recreation and protection of environment or cultural heritage
#•  leisure=nature_reserve - for wildlife or geological / special interest features, managed for conservation, study, research opportunities
#•  boundary=protected_area - includes wilderness, marine protection areas, heritage sites and similar areas
#•  place=square - a town square, typically crossed by streets, also pedestrian areas or rarely green areas
#•  landuse=recreation_ground - an open green space for general recreation (can be public or private)
#•  landuse=village_green - a distinctive part of a village centre
#•  landuse=grass - a smaller area of mown and managed grass
#•  leisure=garden - a place where flowers and other plants are grown
#•  tourism=theme_park - an amusement park where entertainment is provided by rides, game concessions, etc.
#•  tourism=picnic_site

q <- opq(bbox = c(lim$xmin, lim$ymin, lim$xmax, lim$ymax))%>% 
  add_osm_feature(key = "landuse") %>% 
  osmdata_sf()
qosm_sf <- q$osm_polygons
st_crs(qosm_sf) <-4326


q <- opq(bbox = c(lim$xmin, lim$ymin, lim$xmax, lim$ymax))%>% 
add_osm_features (features = list (
        "boundary" = "national_park",
        "boundary" = "protected_area",
        "leisure" = "nature_reserve",
        "place" = "square",
         "landuse" = "recreation_ground",
        "landuse" = "village_green",
         "landuse" = "grass",
        "leisure" = "garden",
        "tourism" = "theme_park",
        "tourism" = "picnic_site"))%>% 
  osmdata_sf()
qosm_sf <- q$osm_polygons
st_crs(qosm_sf) <-4326


vq <- opq(bbox = c(lim$xmin, lim$ymin, lim$xmax, lim$ymax))%>% 
  add_osm_feature(key = "highway") %>% 
  osmdata_sf()
vqosm_sf <- vq$osm_lines 
st_crs(vqosm_sf) <-4326


vqosm_sf <- vqosm_sf|> filter(highway %in% c("primary","secondary","tertiary","residential"))


g <- ggplot() +
      
 geom_raster(
    data = NDVI_df,
    aes(x = x, y = y, fill = ndvi_2023_MX_013006_3)
  ) +
  geom_raster(
    data = NDVI_df2,
    aes(x = x, y = y, fill = ndvi_2023_MX_013007_3)
  ) +
  scale_fill_gradientn(name = "NDVI",colours = cpt(n = 254, "ncl_WhiteGreen"))+
  
  geom_sf(data = qosm_sf, fill= "green",color = "green")+
  geom_sf(data = parquesinegi$geometry, fill= "green",color = "green") +

  geom_sf(data = recortedatasia, fill= "green",color = "green") +
    geom_sf(data = recortedatasiaD, fill= "green",color = "green") +
    geom_sf(data = vqosm_sf, fill = "navy",color="navy",linewidth = 0.05)+
  labs(
    title = "OSM e INEGI SIA parques y jardines en La Laguna",
    subtitle = "01-06-2020",
    caption = "Parques jardines INEGI Osm",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()+
  theme(
    axis.line = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    legend.position = "bottom",
   #legend.position = "none",
        legend.key.width = unit(4, 'cm'), #change legend key width

    legend.title = element_text(size = 11, color = "grey10"),
    legend.text = element_text(size = 10, color = "grey10"),
    panel.grid.major = element_line(color = "white"),
    panel.grid.minor = element_line(color = "white"),
    plot.background = element_rect(fill = "white", color = NA),
    legend.background = element_rect(fill = "white", color = NA),
    panel.border = element_rect(fill = NA, color = "white"))+
  coord_sf(xlim,ylim)

g

#library(ggpubr)
#ggarrange(g,h, common.legend = FALSE)

ggsave("veglocal2.svg", width = 1300, height = 1300, units = "mm", dpi = "retina",limitsize = FALSE)
ggsave("veglocal2.png", width = 300, height = 300, units = "mm", dpi = 300,limitsize = FALSE)