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

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)