1.Introducción

Este es el segundo cuaderno R Notebook en el curso Geomatica Básica 2023 dictada por el ing. Ivan Lizarazo.Aqui ilustraremos cómo hacer mapas temáticos que muestren la participación municipal de los grupos de cultivos más importantes para un el departamento de Huila. Usaremos como fuente principal los archivos csv guardados en el cuaderno EVA, así como un shapefile de municipios obtenidos en clase.

2.Configuración

En este paso, instalaremos y cargaremos las bibliotecas R necesarias.

#ejecute las siguientes líneas desde la línea de comando:
#install.packages('dplyr')
#install.packages('readxl')
#install.packages('sf)
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(readr)

3.Leer archivos relacionados con municipios, cultivos y ciudades.

Antes de comenzar, asegúrese de haber guardado los siguientes archivos en su computadora: 1.Un archivo shapefile con los municipios de su departamento (como se produjo en clase) 2.Un archivo csv con estadísticas EVA 2020 para el grupo de cultivos seleccionado (como se produjo en el primer cuaderno )

Además, es necesario tener un expediente con las principales ciudades de Colombia. Puede descargarlo desde https://simplemaps.com/data/co-cities. Luego, debe descomprimirlo y moverlo a su directorio de trabajo.

Comencemos enumerando los archivos disponibles.

¿Cuáles son los archivos csv guardados en el directorio de trabajo?

list.files("./", pattern=c('csv'))
## [1] "co.csv"                  "huila_frutales_2022.csv"
list.files()
##  [1] "co.csv"                    "ESRI Shapefile.gdb"       
##  [3] "ESRI Shapefile.qmd"        "Geo_Huila.geojson"        
##  [5] "Geo_Huila.qmd"             "H_la_mun.cpg"             
##  [7] "H_la_mun.dbf"              "H_la_mun.prj"             
##  [9] "H_la_mun.qmd"              "H_la_mun.shp"             
## [11] "H_la_mun.shx"              "huila_frutales_2022.csv"  
## [13] "huila_mun_municipios.gpkg" "license.txt"              
## [15] "Mapastematicos.html"       "Mapastematicos.nb.html"   
## [17] "Mapastematicos.Rmd"        "MapaTematico.qgz"         
## [19] "MapaTematico01.qgz"        "MUNICIPIOS.gpkg"          
## [21] "Municipios_de_Huila.cpg"   "Municipios_de_Huila.dbf"  
## [23] "Municipios_de_Huila.prj"   "Municipios_de_Huila.qmd"  
## [25] "Municipios_de_Huila.shp"   "Municipios_de_Huila.shx"  
## [27] "MUNICIPIOS_HUILA.gpkg"

Ahora procedemos a leer los archivos EVA creados en el anterior cuaderno de R

(frutales = read_csv("huila_frutales_2022.csv",show_col_types = FALSE))

Este codigo lee los municipios de el departamento:

(mun.tmp = st_read ("./Municipios_de_Huila.shp"))
## Reading layer `Municipios_de_Huila' from data source 
##   `D:\Sexta Matricula\Geomatica\Clase4-10\Municipios_de_Huila.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 37 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS:  MAGNA-SIRGAS
# comprobar la salida del objeto en el último fragmento y
# cambiar los nombres de los atributos según sus propios datos
mun.tmp %>% select(MPIO_CCDGO, MPIO_CNMBR, AREA) -> municipios
municipios

Ahora leeremos el archivo csv de las ciudades:

#Este archivo se descargó de *simplemaps* como se decribio anteriormente
(cities = read_csv("co.csv"))
## Rows: 82 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): city, country, iso2, admin_name, capital
## dbl (4): lat, lng, population, population_proper
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Tenga en cuenta que hay 82 filas que representan ciudades en Colombia. Tenga en cuenta que el objeto de ciudades es un objeto no espacial (un tibble, en realidad).

Tenga en cuenta también que los valores de las coordenadas se almacenan en los atributos “lng” (latitud) y “lat” (longitud).

Para convertir ciudades en un objeto espacial (un objeto de característica simple, en este caso):

sf.cities <-  st_as_sf(x = cities, 
                        coords = c("lng", "lat"))
sf.cities

¡Tenga en cuenta el CRS que falta!

Agreguemos el CRS usando el código EPSG:

st_crs(sf.cities) <-4326

4. Datos de subconjunto relevantes para el departamento.

Como estamos interesados en un solo departamento, necesitamos crear una unión espacial (spatial join):

sf.cities <- st_transform(sf.cities, st_crs(municipios))
sf.cities.joined <- st_join(sf.cities, municipios, join = st_within)

para mostrar el objeto unido

sf.cities.joined

Tenga en cuenta que obtuvimos un marco de datos sf con cada fila de ciudades adjuntadas con las columnas de nuestros municipios. Las ciudades ubicadas en un departamento diferente tienen muchas NAs.

Ahora filtramos las filas que corresponden a nuestro departamento (en este caso Huila):

huila.cities = dplyr::filter(sf.cities.joined, admin_name=='Huila')
huila.cities

5. Haz un mapa para el grupo de cultivos más importante.

Ahora, haz un mapa coroplético de estos datos. Usaremos las bibliotecas tmap, ggplot2 y classInt.

# ejecuta las siguientes líneas desde la ventana de comandos:
#install.packages("tmap")
#install.packages("ggplot2")
#
library(ggplot2)
library(terra)
## terra 1.7.46
library(tmap)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(classInt)

Recordemos que el objeto municipios no tiene atributos EVA. En cambio, el objeto frutales, que contiene estadísticas de cultivos, es un objeto no espacial. Nuestra siguiente tarea es unir el objeto de estadísticas a los objetos espaciales para tener los datos relevantes en un solo objeto.

Para poder realizar la unión, necesitamos una clave compartida, es decir, un atributo común. En este caso lo tenemos en ambos objetos, ese es el código municipal. Sin embargo, tenga en cuenta que tanto sus nombres como sus tipos de datos son diferentes.

class(frutales$COD_MUN)
## [1] "numeric"
class(municipios$MPIO_CCDGO)
## [1] "character"

Por lo tanto, necesitamos arreglar esto:

frutales$COD_MUN = as.character(frutales$COD_MUN)

Ahora estamos listos para unirnos. Es posible que necesites revisar tus notas de clase para entender lo que sucede aquí:

munic_frutales = left_join(municipios, frutales, by = c("MPIO_CCDGO" = "COD_MUN"))

Cuales son los resultados ?

munic_frutales

Tenga en cuenta que munic_frutales incluye ahora el atributo max_prod que representa la cantidad total de toneladas producidas por cultivos frutales en 2020.

El siguiente código personaliza los colores de los municipios. Más información aqui (https://stackoverflow.com/questions/57177312/trying-to-plot-in-tmap-shapefile-with-attribute).