Este cuaderno muestra como agregar atributos a un objeto simple espacial (simple feature).
primero, cargamos la librería sf
library(sf)
Luego, vamos a leer un shape file de los municipios del departamento de Bolívar:
(BOL_mun <- read_sf("mbolivar.shp"))
Simple feature collection with 36 features and 6 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -75.79764 ymin: 6.979202 xmax: -73.769 ymax: 10.80153
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
Ahora vamos a identificar la clase de esa variable (BOL_mun)
class(BOL_mun)
[1] "sf" "tbl_df" "tbl" "data.frame"
Luego vamos a identificar los atributos de “BOL_mun”
attributes(BOL_mun)
$names
[1] "ISO" "NAME_0" "NAME_1" "NAME_2" "TYPE_2" "ENGTYPE_2" "geometry"
$row.names
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
$class
[1] "sf" "tbl_df" "tbl" "data.frame"
$sf_column
[1] "geometry"
$agr
ISO NAME_0 NAME_1 NAME_2 TYPE_2 ENGTYPE_2
<NA> <NA> <NA> <NA> <NA> <NA>
Levels: constant aggregate identity
Ahora vamos a conocer el sistema de referencia de coordenadas (CRS)
st_crs(BOL_mun)
Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"
Retrieving bounding box
st_bbox(BOL_mun)
xmin ymin xmax ymax
-75.797638 6.979202 -73.768997 10.801529
Ahora vamos a instalar algunos paquetes adicionales, acordarse que el nombre de las librerias va en comillas
###install.packages("lwgeom")
###install.packages("units")
library(lwgeom)
library(units)
(BOL_mun$area <- st_area(BOL_mun))
Units: [m^2]
[1] 1518018318 271625128 542778128 406004232 365903883 1085708052 758102169 1363977282 497471890
[10] 811980212 352170859 147956277 465248779 417747848 461708919 603146279 2025289041 3121393378
[19] 1618697738 403204500 191907613 359660556 438286337 640093298 387286695 1141735828 175030538
[28] 128840937 2488265513 1782555699 120915336 127105232 182761350 183267145 157324656 291402214
Cambio de unidades del área, ya que la anterior estaba en m2, ahora la vamos a pasar a km2
set_units(BOL_mun$area, km^2)
Units: [km^2]
[1] 1518.0183 271.6251 542.7781 406.0042 365.9039 1085.7081 758.1022 1363.9773 497.4719 811.9802
[11] 352.1709 147.9563 465.2488 417.7478 461.7089 603.1463 2025.2890 3121.3934 1618.6977 403.2045
[21] 191.9076 359.6606 438.2863 640.0933 387.2867 1141.7358 175.0305 128.8409 2488.2655 1782.5557
[31] 120.9153 127.1052 182.7614 183.2671 157.3247 291.4022
Acontinuación vamos a usar el Censo agropecuario
library(readxl)
#xls filies
(my_data <- read_excel("CENSO.xls"))
En esta parte se va a escoger mi departamento Bolívar, donde se muestra una tabla con los atributos (datos de censo para nuestro departamento) ### Joining censodata to simple feature object:
(BOLIVAR_DATA <- my_data %>% filter(DEPTO == "Bolívar"))
(BOLIVAR_NUEVO = rename(BOLIVAR_DATA, NAME_2 = MUNICIPIO))
Ahora es tiempo de hacer una unión “left_join”:
(bolivar_censo <- left_join(BOL_mun, BOLIVAR_NUEVO))
Joining, by = "NAME_2"
Simple feature collection with 36 features and 14 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -75.79764 ymin: 6.979202 xmax: -73.769 ymax: 10.80153
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
De los datos de hacer el “left_join”, es el mismo código del cuaderno 1, solo cambia la primera vaiable censo
bolivar_censo %>% mapview(zcol = "NAME_2", legend = TRUE, col.regions = sf.colors)
Ahora vamos a ordenar las municipalidades basados en el censo. Primero por el área de bosque en hectáreas… Se organiza en forma descendente, (NA es que no tiene datos - nueva variable)
(bolivar_clean <- bolivar_censo %>% arrange(desc(HA_BOSQUE)) %>% select(NAME_2, MCOD, HA_BOSQUE, HA_AGROP, HA_NO_AGROP, HA_USOS))
Simple feature collection with 36 features and 6 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -75.79764 ymin: 6.979202 xmax: -73.769 ymax: 10.80153
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs