análisis_Espacial_04

Lineas de visión

determina la visibilidad de las lineas de visión sobre obstáculos formados por una superfiie y un data set opcional, se genera ocn la extensión 3D

análisis de cambio de uso de suelo en colima a nivel estatal y por cuenca

Clases pára reclasificar: permanencia natural, permanencia antrópica, deforestación, recuepración, falso cambio.

Productos a elaborar:

• Tabla de cambios con superficies en hectareas estatal y por cuenca
• Mapa de salida con sombreado del MDT
• Limite de cuencas

library(tidyverse)

ent = sf::st_read("C:/Users/DELL/Desktop/Marco_Geo_2020/889463807469_s/MG_2020_Integrado/conjunto_de_datos/00ent.shp",quiet=T) %>% 
  janitor::clean_names() %>% 
  filter(cve_ent=="06") %>% 
  sf::st_transform(crs=32613)

bbox_colima = sf::st_read("../01_input/bbox_colima.gpkg",quiet=T)

ent = sf::st_intersection(ent %>% sf::st_transform(crs=4326),bbox_colima)



serie_5 = sf::st_read("../01_input/análisis espacial/usv250s5ugw/usv250s5ugw.shp",quiet=T) %>%   sf::st_transform(crs=4326) %>% janitor::clean_names()

serie_7 = sf::st_read("../01_input/análisis espacial/usv250s7gw/usv250s7gw.shp",quiet=T) %>%   sf::st_transform(crs=4326) %>% janitor::clean_names()

Corte

sf::sf_use_s2(F)
serie_5 = serie_5 %>% 
  sf::st_join(ent %>% select(1),join = sf::st_intersects) %>% 
  filter(!is.na(cvegeo)) %>% 
  sf::st_intersection(ent %>% select(1))

serie_7 = serie_7 %>% 
  sf::st_join(ent %>% select(1),join = sf::st_intersects) %>% 
  filter(!is.na(cvegeo)) %>% 
  sf::st_intersection(ent %>% select(1))

Visualización

library(leaflet)

leaflet() %>% 
  addProviderTiles(providers$CartoDB) %>% 
  addPolygons(data = ent,weight = 2,fill = NA,color = "black",opacity = 1) 

Get base map

mapa_base = ggplot()+
  basemapR::base_map(bbox = sf::st_bbox(ent),basemap = "voyager",nolabels = T,increase_zoom = 5)

mapa_5 = mapa_base+
  geom_sf(data=serie_5,aes(fill=descripcio),show.legend = F)+
  geom_sf(data=ent,fill=NA,color="black")+
  labs(title = "Serie 5")+
  ggthemes::theme_map()

mapa_7 = mapa_base+
  geom_sf(data=serie_7,aes(fill=descripcio),show.legend = F)+
  geom_sf(data=ent,fill=NA,color="black")+
  labs(title = "Serie 7")+
  ggthemes::theme_map()

cowplot::plot_grid(mapa_5,mapa_7)

Reclasificación

serie_5 = serie_5 %>% 
  mutate(descripcio=tolower(descripcio)) %>% 
  mutate(reclas=case_when(
    str_detect(descripcio,"bosque")~"Bosque",
    str_detect(descripcio,"agricultura")~"Cultivo",
    str_detect(descripcio,"secundaria")~"Bosque",
    str_detect(descripcio,"selva")~"Selva",
    str_detect(descripcio,"cuerpo de agua")~"Otras coberturas",
    str_detect(descripcio,"zona urbana")~"Otras coberturas",
    str_detect(descripcio,"asentamientos humanos")~"Otras coberturas",
    str_detect(descripcio,"manglar")~"Vegetación hidrófila",
    str_detect(descripcio,"tular")~"Vegetación hidrófila",
    str_detect(descripcio,"vegetación halófila hidrófila")~"Vegetación hidrófila",
    str_detect(descripcio,"palmar")~"Otros tipos de vegetación",
    str_detect(descripcio,"pastizal")~"Pastizal",
    str_detect(descripcio,"sabana")~"Pastizal",
    str_detect(descripcio,"pradera")~"Pastizal",
    T~"Otros tipos de vegetación"
  ))

serie_7= serie_7 %>% 
  mutate(descripcio=tolower(descripcio)) %>% 
  mutate(reclas=case_when(
    str_detect(descripcio,"bosque")~"Bosque",
    str_detect(descripcio,"agricultura")~"Cultivo",
    str_detect(descripcio,"secundaria")~"Bosque",
    str_detect(descripcio,"selva")~"Selva",
    str_detect(descripcio,"cuerpo de agua")~"Otras coberturas",
    str_detect(descripcio,"zona urbana")~"Otras coberturas",
    str_detect(descripcio,"asentamientos humanos")~"Otras coberturas",
    str_detect(descripcio,"manglar")~"Vegetación hidrófila",
    str_detect(descripcio,"tular")~"Vegetación hidrófila",
    str_detect(descripcio,"vegetación halófila hidrófila")~"Vegetación hidrófila",
    str_detect(descripcio,"palmar")~"Otros tipos de vegetación",
    str_detect(descripcio,"pastizal")~"Pastizal",
    str_detect(descripcio,"sabana")~"Pastizal",
    str_detect(descripcio,"pradera")~"Pastizal",
    T~"Otros tipos de vegetación"
  ))

area_5 = serie_5 %>% 
  sf::st_transform(crs=32614) %>% 
  sf::st_make_valid() %>% 
  group_by(reclas) %>% 
  summarise(total = units::set_units(sum(sf::st_area(geometry)), "km^2"))

area_7 = serie_7 %>% 
  sf::st_transform(crs=32614) %>% 
  sf::st_make_valid() %>% 
  group_by(reclas) %>% 
  summarise(total = units::set_units(sum(sf::st_area(geometry)), "km^2"))

area_5 %>% 
  as_tibble() %>% 
  select(reclas,total_5=total) %>% 
  left_join(area_7 %>% select(reclas,total_7=total) %>% as_tibble(),by="reclas") %>% 
  pivot_longer(cols = c(total_5,total_7),names_to = "serie",values_to = "area_km2") %>% 
  mutate(area_km2=as.numeric(area_km2)) %>% 
  ggplot(aes(reclas,area_km2,fill=serie))+
  geom_col(position = "dodge2")+
  facet_wrap(~reclas,scales = "free")+
  theme_bw()+
  theme(legend.position = "top")+
  scale_fill_manual(values = wesanderson::wes_palette("Darjeeling1"))

Ahora se verá el cambio de suelo

Vamos a hacer una intersección las dos geometrías por lo que tendremos en una misma capa, o sea cruzar, el resultado tendrá dos tablas de atributos previas. había intentado con sf, pero es muy lento por lo que usaré terra :)

terra_5 = serie_5 %>%
  rename(reclas_serie_5=reclas) %>% 
  select(reclas_serie_5)

terra_7 = serie_7 %>%
  rename(reclas_serie_7=reclas) %>% 
  select(reclas_serie_7)

intersección con terra::

# interseccion= terra::intersect(terra_5,terra_7)

llamado proceso qgis

siento que se tarda mucho, entonces creo que con qgis podría ser más rapido, por lo que llamaré el argumento y los elementos para hacer la intersección. con la función ger argument obtenemos los argumentos para llamar al proceso

qgisprocess::qgis_get_argument_specs("native:intersection")

## Attempting to load the package cache ... Success!

## # A tibble: 7 × 6
##   name    description qgis_type default_value available_values acceptable_values
##   <chr>   <chr>       <chr>     <list>        <list>           <list>           
## 1 INPUT   Input layer source    <NULL>        <NULL>           <chr [1]>        
## 2 OVERLAY Overlay la… source    <NULL>        <NULL>           <chr [1]>        
## 3 INPUT_… Input fiel… field     <NULL>        <NULL>           <chr [2]>        
## 4 OVERLA… Overlay fi… field     <NULL>        <NULL>           <chr [2]>        
## 5 OVERLA… Overlay fi… string    <NULL>        <NULL>           <chr [3]>        
## 6 OUTPUT  Intersecti… sink      <NULL>        <NULL>           <chr [1]>        
## 7 GRID_S… Grid size   number    <NULL>        <NULL>           <chr [3]>

proceso_qgis = qgisprocess::qgis_run_algorithm(
      algorithm = "native:intersection",
      INPUT =terra_5,
      OVERLAY=terra_7
      )

## Argument `INPUT_FIELDS` is unspecified (using QGIS default value).

## Argument `OVERLAY_FIELDS` is unspecified (using QGIS default value).

## Argument `OVERLAY_FIELDS_PREFIX` is unspecified (using QGIS default value).

## Using `OUTPUT = qgis_tmp_vector()`

## Argument `GRID_SIZE` is unspecified (using QGIS default value).

resultado = sf::st_read(proceso_qgis$OUTPUT)

## Reading layer `file40bc5a224ba7' from data source 
##   `C:\Users\DELL\AppData\Local\Temp\RtmpmM74bk\file40bc7c822370\file40bc5a224ba7.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 3381 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -104.6906 ymin: 18.68415 xmax: -103.4863 ymax: 19.51252
## Geodetic CRS:  WGS 84

clasificacion = resultado %>% 
  as_tibble() %>% 
  count(reclas_serie_5,reclas_serie_7,sort = T)

clasificacion

## # A tibble: 45 × 3
##    reclas_serie_5   reclas_serie_7       n
##    <chr>            <chr>            <int>
##  1 Bosque           Bosque             600
##  2 Cultivo          Cultivo            390
##  3 Cultivo          Bosque             357
##  4 Bosque           Pastizal           261
##  5 Bosque           Cultivo            214
##  6 Pastizal         Bosque             176
##  7 Cultivo          Otras coberturas   129
##  8 Otras coberturas Cultivo            122
##  9 Cultivo          Pastizal           117
## 10 Pastizal         Pastizal            99
## # ℹ 35 more rows

clasificacion = readxl::read_excel("../01_input/clasificacion_llm.xlsx")

clasificacion = clasificacion %>% 
  mutate(llave=paste0(serie_5,serie_7))

clasificacion

## # A tibble: 45 × 4
##    serie_5          serie_7          clas_cambio           llave                
##    <chr>            <chr>            <chr>                 <chr>                
##  1 Bosque           Bosque           Permanencia espacial  BosqueBosque         
##  2 Cultivo          Cultivo          Permanencia espacial  CultivoCultivo       
##  3 Cultivo          Bosque           Recuperación          CultivoBosque        
##  4 Bosque           Pastizal         Deforestación         BosquePastizal       
##  5 Bosque           Cultivo          Deforestación         BosqueCultivo        
##  6 Pastizal         Bosque           Recuperación          PastizalBosque       
##  7 Cultivo          Otras coberturas Permanencia antrópica CultivoOtras cobertu…
##  8 Otras coberturas Cultivo          Permanencia antrópica Otras coberturasCult…
##  9 Cultivo          Pastizal         Permanencia antrópica CultivoPastizal      
## 10 Pastizal         Pastizal         Permanencia espacial  PastizalPastizal     
## # ℹ 35 more rows

resultado_final = resultado %>% 
    mutate(llave=paste0(reclas_serie_5,reclas_serie_7)) %>% 
  left_join(clasificacion %>% select(llave,clas_cambio),by="llave")

mapa_base+
  geom_sf(data=resultado_final,aes(fill=clas_cambio,col=clas_cambio),alpha=.7)

mapa_base+
  geom_sf(data=resultado_final,aes(fill=clas_cambio,col=clas_cambio),alpha=.7)+
  facet_wrap(~clas_cambio)+
  theme(legend.position = "top")

resultado_final$clas_cambio=factor(resultado_final$clas_cambio)

color_piso = colorFactor(palette = "viridis",domain = resultado_final$clas_cambio)

table(resultado_final$clas_cambio)

## 
##         Deforestación          Falso cambio Permanencia antrópica 
##                   674                   201                   493 
##  Permanencia espacial          Recuperación 
##                  1257                   756

leaflet() %>% 
  addProviderTiles(providers$Esri.WorldImagery,group = "Satelite") %>% 
  addProviderTiles(providers$OpenStreetMap,group = "OSM") %>% 
  addPolygons(data=resultado_final,
              color = ~color_piso(clas_cambio),
              fillColor = ~color_piso(clas_cambio),
              group = ~clas_cambio,
              weight = 1,
              fillOpacity = .7) %>% 
  addLayersControl(
    baseGroups = c("Satelite","OSM"),
    overlayGroups = c(unique(resultado_final$clas_cambio)),
    options = layersControlOptions(collapsed=F)
    )

resultado_final = resultado_final %>% 
  mutate(label=paste0("Serie 5: ",reclas_serie_5,"| Serie 7:", reclas_serie_7,"| Cambio:", clas_cambio))

leaflet() %>% 
  addProviderTiles(providers$Esri.WorldImagery, group = "Satelite") %>% 
  addProviderTiles(providers$OpenStreetMap, group = "OSM") %>% 
  addPolygons(data = resultado_final,
              color = ~color_piso(clas_cambio),
              fillColor = ~color_piso(clas_cambio),
              group = ~clas_cambio,
              weight = 1,
              fillOpacity = .7,
              highlightOptions = highlightOptions(
                weight = 2,
                color = "#FF0000",
                fillOpacity = 0.7,
                bringToFront = TRUE
              ),
              label = ~label) %>% 
  addLayersControl(
    baseGroups = c("Satelite", "OSM"),
    overlayGroups = c(unique(resultado_final$clas_cambio)),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
  addLegend(pal = color_piso,
            values = resultado_final$clas_cambio,
            opacity = 1,
            title = "Cambios de uso de suelo"
            )

CEM COLIMA

• Leemos el archivo
• reproyectamos a zona 13
• visualizamos

cem = terra::rast("../01_input/CEM_V3_20170619_R15_E06_TIF/Colima_r15m.tif")

utm_zona14n = "+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs"

cem= terra::project(cem, utm_zona14n)

terra::plet(cem)

orientacion <- terra::terrain(x = cem, v = "aspect")
terra::plet(orientacion)

Hilshade

hillshade <- terra::shade(
  slope = terra::terrain(cem, "slope"),
  aspect = terra::terrain(cem, "aspect"),
  angle = 45,  # ángulo de elevación del sol
  direction = 315,  # dirección de la fuente de luz (noroeste)  
  normalize = TRUE
)

terra::plet(hillshade)

Lineas de visión se hará con Qgis

qgisprocess::qgis_algorithms() %>% 
  filter(grepl("sight|visibility|viewshed", algorithm))

## # A tibble: 4 × 24
##   provider provider_title algorithm                 algorithm_id algorithm_title
##   <chr>    <chr>          <chr>                     <chr>        <chr>          
## 1 gdal     GDAL           gdal:viewshed             viewshed     Viewshed       
## 2 grass    GRASS          grass:r.viewshed          r.viewshed   r.viewshed     
## 3 grass    GRASS          grass:v.net.visibility    v.net.visib… v.net.visibili…
## 4 sagang   SAGA Next Gen  sagang:visibilityanalysis visibilitya… Visibility ana…
## # ℹ 19 more variables: provider_can_be_activated <lgl>,
## #   provider_is_active <lgl>, provider_long_name <chr>, provider_version <chr>,
## #   provider_warning <chr>, can_cancel <lgl>, deprecated <lgl>, group <chr>,
## #   has_known_issues <lgl>, help_url <chr>, requires_matching_crs <lgl>,
## #   short_description <chr>, tags <list>, default_raster_file_extension <chr>,
## #   default_vector_file_extension <chr>,
## #   supported_output_raster_extensions <list>, …

al identificar el algoritmo ahora buscaremos sus parámetros

qgisprocess::qgis_get_argument_specs(
  "grass:r.viewshed"
  )

## # A tibble: 16 × 6
##    name   description qgis_type default_value available_values acceptable_values
##    <chr>  <chr>       <chr>     <list>        <list>           <list>           
##  1 input  Elevation   raster    <NULL>        <NULL>           <chr [1]>        
##  2 coord… Coordinate… point     <chr [1]>     <NULL>           <chr [1]>        
##  3 obser… Viewing el… number    <dbl [1]>     <NULL>           <chr [3]>        
##  4 targe… Offset for… number    <dbl [1]>     <NULL>           <chr [3]>        
##  5 max_d… Maximum vi… number    <dbl [1]>     <NULL>           <chr [3]>        
##  6 refra… Refraction… number    <dbl [1]>     <NULL>           <chr [3]>        
##  7 memory Amount of … number    <dbl [1]>     <NULL>           <chr [3]>        
##  8 -c     Consider e… boolean   <lgl [1]>     <NULL>           <chr [4]>        
##  9 -r     Consider t… boolean   <lgl [1]>     <NULL>           <chr [4]>        
## 10 -b     Output for… boolean   <lgl [1]>     <NULL>           <chr [4]>        
## 11 -e     Output for… boolean   <lgl [1]>     <NULL>           <chr [4]>        
## 12 output Intervisib… rasterDe… <NULL>        <NULL>           <chr [1]>        
## 13 GRASS… GRASS GIS … extent    <NULL>        <NULL>           <chr [2]>        
## 14 GRASS… GRASS GIS … number    <dbl [1]>     <NULL>           <chr [3]>        
## 15 GRASS… Output Ras… string    <NULL>        <NULL>           <chr [3]>        
## 16 GRASS… Output Ras… string    <NULL>        <NULL>           <chr [3]>

Ahora el modelo de view

puntos_colima = sf::st_read("../01_input/puntos_vision_colima.gpkg",quiet=T)
puntos_colima=puntos_colima %>% sf::st_transform(crs=32613)