#install.packages("sf")
#install.packages("tmap")
#install.packages("ggthemes")
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(tidyr)
library(tmap)
## Warning: package 'tmap' was built under R version 4.3.3
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(sf)
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(ggplot2)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.3.3
countries <- st_read("C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/World_Countries/World_Countries.shp")
## Reading layer `World_Countries' from data source 
##   `C:\Users\user\OneDrive - Softkey Digital\Aldo Mendoza Cotera\2024-1\Herramientas cuantitativas\Tarea 2\introgeodf\maps\World_Countries\World_Countries.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 252 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6236
## Geodetic CRS:  WGS 84

Introducción

Tipo de la data

class(countries)
## [1] "sf"         "data.frame"

Dimensión de la data

dim(countries)
## [1] 252   2

Nombres

names(countries)
## [1] "COUNTRY"  "geometry"

Contexto

head(as.data.frame(countries))

Valores vacíos

countries[apply(is.na(countries), 1, any), ]
## Warning in `[<-.data.frame`(`*tmp*`, is_list, value = list(`2` =
## "<s_GEOMET>")): replacement element 1 has 1 row to replace 0 rows

Tipo

str(countries)
## Classes 'sf' and 'data.frame':   252 obs. of  2 variables:
##  $ COUNTRY : chr  "Aruba (Netherlands)" "Antigua and Barbuda" "Afghanistan" "Algeria" ...
##  $ geometry:sfc_MULTIPOLYGON of length 252; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:11, 1:2] -69.9 -69.9 -70.1 -70.1 -70 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
##   ..- attr(*, "names")= chr "COUNTRY"

Mapa del mundo sin ríos y ciudades

# Trazar el mapa utilizando ggplot2
ggplot(data = countries) +
  geom_sf(fill = "azure", color = "black", size = 0.1) +
  labs(x = "Longitud", y = "Latitud")

Ejercicio 1

Archivos de ríos y ciudades

rivers <- st_read("C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/World_Hydrography/World_Hydrography.shp")
## Reading layer `World_Hydrography' from data source 
##   `C:\Users\user\OneDrive - Softkey Digital\Aldo Mendoza Cotera\2024-1\Herramientas cuantitativas\Tarea 2\introgeodf\maps\World_Hydrography\World_Hydrography.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 2 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -164.6246 ymin: -36.90435 xmax: 161.0658 ymax: 72.23418
## Geodetic CRS:  WGS 84
cities <- st_read("C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/World_Cities/World_Cities.shp")
## Reading layer `World_Cities' from data source 
##   `C:\Users\user\OneDrive - Softkey Digital\Aldo Mendoza Cotera\2024-1\Herramientas cuantitativas\Tarea 2\introgeodf\maps\World_Cities\World_Cities.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 610 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -165.27 ymin: -53.15 xmax: 177.1302 ymax: 78.2
## Geodetic CRS:  WGS 84

Mapa de ríos

# Trazar los ríos utilizando ggplot2
ggplot(data = rivers) +
  geom_sf(color = "blue", size = 1, linetype = "dotted")

Mapa de ciudades

# Trazar las ciudades utilizando ggplot2
ggplot(data = cities) +
  geom_sf(color = "red", size = 1, alpha = 0.3, shape = 19)

Comprobar si los CRS son los mismos

# Comprobando si los CRS son iguales
if (st_crs(countries) == st_crs(cities)) {
  print("Los CRS de countries y cities son iguales.")
} else {
  print("Los CRS de countries y cities son diferentes.")
}
## [1] "Los CRS de countries y cities son iguales."

Mapa del mundo

# Trazar el mapa base (countries)
map_base <- ggplot(data = countries) +
  geom_sf(fill = "grey90", color = "black", size = 0.1) +
  theme_light() +
  theme(aspect.ratio = 1)

# Trazar ríos sobre el mapa base
map_with_rivers <- map_base +
  geom_sf(data = rivers, color = "blue", size = 0.4)

# Trazar ciudades sobre el mapa con ríos
final_map <- map_with_rivers +
  geom_sf(data = cities, color = "red", size = 1, alpha = 0.7, shape = 19)

final_map

Ejercicio 2

Escribir datos de países, ríos y ciudades

# Escribir countries en el archivo worldMap.gpkg, capa 'countryBorders'
st_write(countries, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "countryBorders", driver = "GPKG")
## Writing layer `countryBorders' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 252 features with 1 fields and geometry type Multi Polygon.
# Escribir rivers en el archivo worldMap.gpkg, capa 'riverLines'
st_write(rivers, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "riverLines", driver = "GPKG")
## Writing layer `riverLines' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 98 features with 2 fields and geometry type Multi Line String.
# Escribir cities en el archivo worldMap.gpkg, capa 'cityPoints'
st_write(cities, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "cityPoints", driver = "GPKG")
## Writing layer `cityPoints' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 610 features with 3 fields and geometry type Point.
# Filtrar las filas donde el país es 'Peru'

peru <- filter(countries, COUNTRY == 'Peru')

Recortando para Peru

# Realizar el recorte de las ciudades
citiesPeru_clipped <- st_intersection(cities, peru)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# Realizar el recorte de los ríos
riversPeru_clipped <- st_intersection(rivers, peru)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

Mapa Peru

# Tamaño de la figura
plot_size <- 5

# Trazar el mapa base de Brasil
base <- ggplot(data = peru) +
  geom_sf(fill = "grey90", color = "black", size = 0.4) +
  theme_light() +
  theme(legend.position = "none") +
  coord_sf()

# Trazar ciudades sobre el mapa base de Brasil
base + 
  geom_sf(data = citiesPeru_clipped, color = "red", size = 2) +
  geom_sf(data = riversPeru_clipped, color = "blue", size = 1)
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Ejercicio 3

Proyeccion de mapa

st_crs(peru)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
#centroide
st_centroid(peru)
## Warning: st_centroid assumes attributes are constant over geometries
# Transformar Brasil al nuevo CRS (5641)
st_transform(peru, 5641) %>%
  st_crs()
## Coordinate Reference System:
##   User input: EPSG:5641 
##   wkt:
## PROJCRS["SIRGAS 2000 / Brazil Mercator",
##     BASEGEOGCRS["SIRGAS 2000",
##         DATUM["Sistema de Referencia Geocentrico para las AmericaS 2000",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4674]],
##     CONVERSION["Petrobras Mercator",
##         METHOD["Mercator (variant B)",
##             ID["EPSG",9805]],
##         PARAMETER["Latitude of 1st standard parallel",-2,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Longitude of natural origin",-43,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",5000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",10000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Oil and gas exploration and production."],
##         AREA["Brazil - offshore - equatorial margin."],
##         BBOX[-5.74,-51.64,7.04,-32.43]],
##     ID["EPSG",5641]]
# Calcular el centroide de Peru
st_transform(peru, 5641) %>%
  st_centroid()
## Warning: st_centroid assumes attributes are constant over geometries
# Trazar el mapa de Peru con ggplot2 y los centroides de los polígonos
peru %>%
  st_transform(5641) %>%
  ggplot() +
  geom_sf() +
  geom_sf(data = . %>% st_centroid(), color = "red")
## Warning: st_centroid assumes attributes are constant over geometries

# Mantengamos la versión proyectada para todos nuestros mapas:

peru_5641 <- st_transform(peru, 5641)

# Transformar ciudades recortadas al CRS de Brasil transformado

cities_peru_5641 <- st_transform(citiesPeru_clipped, st_crs(peru_5641))

# Transformar ríos recortados al CRS de Brasil transformado
rivers_peru_5641 <- st_transform(riversPeru_clipped, st_crs(peru_5641))
# Guardando 

# Escribir las capas en el archivo GeoPackage

st_write(peru_5641, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "country", driver = "GPKG")
## Writing layer `country' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 1 features with 1 fields and geometry type Multi Polygon.
st_write(cities_peru_5641, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "cities", driver = "GPKG")
## Writing layer `cities' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 8 features with 4 fields and geometry type Point.
st_write(rivers_peru_5641, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "rivers", driver = "GPKG")
## Writing layer `rivers' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 5 features with 3 fields and geometry type Unknown (any).
peru_5641_centroid <- st_centroid(peru_5641)
## Warning: st_centroid assumes attributes are constant over geometries
st_write(peru_5641_centroid, "C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg", layer = "centroid", driver = "GPKG")
## Writing layer `centroid' to data source 
##   `C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/introgeodf/maps/worldMap.gpkg' using driver `GPKG'
## Writing 1 features with 1 fields and geometry type Point.

Creando data espacial

Seguramente obtendrá líneas y polígonos como mapas, pero puede que ese no sea el caso con los puntos. Déjame descargar un archivo CSV con información de los aeropuertos de Brasil desde este sitio web, lo guardaré en mi carpeta de datos:

# Cargar el paquete readr si aún no está cargado

# install.packages("readr")

library(readr)

# Leer el archivo CSV
infoairports <- read_csv("C:/Users/user/OneDrive - Softkey Digital/Aldo Mendoza Cotera/2024-1/Herramientas cuantitativas/Tarea 2/pe-airports.csv")
## Rows: 203 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (23): id, ident, type, name, latitude_deg, longitude_deg, elevation_ft, ...
## 
## ℹ 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.
# Seleccionar las filas específicas
infoairports[c(1:4, nrow(infoairports) - 3:nrow(infoairports)), ]

Esto necesita algo de limpieza:

# Primero la primera fila

# Omitir la primera fila
infoairports <- infoairports[-1, ]

# Restablecer el índice
row.names(infoairports) <- NULL

# Imprimir las primeras filas
head(infoairports)

Mantenemos las columnas necesarias

# Obtener los nombres de las columnas
names(infoairports)
##  [1] "id"                "ident"             "type"             
##  [4] "name"              "latitude_deg"      "longitude_deg"    
##  [7] "elevation_ft"      "continent"         "country_name"     
## [10] "iso_country"       "region_name"       "iso_region"       
## [13] "local_region"      "municipality"      "scheduled_service"
## [16] "gps_code"          "iata_code"         "local_code"       
## [19] "home_link"         "wikipedia_link"    "keywords"         
## [22] "score"             "last_updated"
# Seleccionar columnas específicas
keep <- c('name', 'type', 'latitude_deg', 'longitude_deg', 'elevation_ft', 'region_name', 'municipality')

infoairports <- infoairports[, keep]

# Imprimir el marco de datos resultante
print(infoairports)
## # A tibble: 202 × 7
##    name   type  latitude_deg longitude_deg elevation_ft region_name municipality
##    <chr>  <chr> <chr>        <chr>         <chr>        <chr>       <chr>       
##  1 Jorge… larg… -12.0219     -77.114305    113          Lima Region Lima        
##  2 Aleja… larg… -13.5356998… -71.93879699… 10860        Cuzco Regi… Cusco       
##  3 Rodrí… medi… -16.3411006… -71.58309936… 8405         Arequipa R… Arequipa    
##  4 Inca … medi… -15.4671001… -70.158203125 12552        Puno Region Juliaca     
##  5 Maria… medi… -14.854      -74.961502    1860         Ica Region  Nazca       
##  6 Coron… medi… -3.78474     -73.3088      306          Loreto Reg… Iquitos     
##  7 Capit… medi… -8.08141040… -79.10880279… 106          La Liberta… Trujillo    
##  8 Air F… medi… -6.78748     -79.828102    97           Lambayeque… Chiclayo    
##  9 Capit… medi… -5.20574998… -80.61640167… 120          Piura Regi… Piura       
## 10 Coron… medi… -18.0533008… -70.27580261… 1538         Tacna Regi… Tacna       
## # ℹ 192 more rows
# Obtener información sobre el marco de datos
str(infoairports)
## tibble [202 × 7] (S3: tbl_df/tbl/data.frame)
##  $ name         : chr [1:202] "Jorge Chávez International Airport" "Alejandro Velasco Astete International Airport" "Rodríguez Ballón International Airport" "Inca Manco Capac International Airport" ...
##  $ type         : chr [1:202] "large_airport" "large_airport" "medium_airport" "medium_airport" ...
##  $ latitude_deg : chr [1:202] "-12.0219" "-13.535699844400002" "-16.3411006927" "-15.467100143432617" ...
##  $ longitude_deg: chr [1:202] "-77.114305" "-71.9387969971" "-71.5830993652" "-70.158203125" ...
##  $ elevation_ft : chr [1:202] "113" "10860" "8405" "12552" ...
##  $ region_name  : chr [1:202] "Lima Region" "Cuzco Region" "Arequipa Region" "Puno Region" ...
##  $ municipality : chr [1:202] "Lima" "Cusco" "Arequipa" "Juliaca" ...

Formateamos

# Convertir columnas específicas a numérico
numericCols <- c('latitude_deg', 'longitude_deg', 'elevation_ft')

infoairports <- infoairports %>%
  mutate_at(vars(numericCols), as.numeric)
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(numericCols)
## 
##   # Now:
##   data %>% select(all_of(numericCols))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Obtener información sobre el marco de datos
str(infoairports)
## tibble [202 × 7] (S3: tbl_df/tbl/data.frame)
##  $ name         : chr [1:202] "Jorge Chávez International Airport" "Alejandro Velasco Astete International Airport" "Rodríguez Ballón International Airport" "Inca Manco Capac International Airport" ...
##  $ type         : chr [1:202] "large_airport" "large_airport" "medium_airport" "medium_airport" ...
##  $ latitude_deg : num [1:202] -12 -13.5 -16.3 -15.5 -14.9 ...
##  $ longitude_deg: num [1:202] -77.1 -71.9 -71.6 -70.2 -75 ...
##  $ elevation_ft : num [1:202] 113 10860 8405 12552 1860 ...
##  $ region_name  : chr [1:202] "Lima Region" "Cuzco Region" "Arequipa Region" "Puno Region" ...
##  $ municipality : chr [1:202] "Lima" "Cusco" "Arequipa" "Juliaca" ...
# Graficar los aeropuertos
ggplot(infoairports, aes(x = longitude_deg, y = latitude_deg)) +
  geom_point() +
  coord_equal() +  # Para asegurar que las escalas de los ejes sean iguales
  theme_minimal()  # Para un aspecto limpio del gráfico

### Peru: rios, ciudades y aeropuertos

# Vamos a graficar

# Agregar los aeropuertos al gráfico
map_with_airports <- base +
  geom_sf(data = citiesPeru_clipped, color = "red", size = 2) +
  geom_sf(data = riversPeru_clipped, color = "blue", size = 1) +
  geom_point(data = infoairports, aes(x = longitude_deg, y = latitude_deg), color = "black", size = 2)
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
# Visualizar el gráfico
print(map_with_airports)