Este cuaderno nos muestra como agregar atributos a un simple feature espacial.
Primero cargamos la libreria sf
library(sf)
library(tidyverse)Luego, leemos un shapefile de municipios en los departamentos de Quindio, Risaralda y Caldas (Eje cafetero).
(Eje_cafetero_mun<-read_sf("Eje_cafetero_mun.shp"))Simple feature collection with 51 features and 6 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defsAhora se identifica que tipo de clase es
class(Eje_cafetero_mun)[1] "sf" "tbl_df" "tbl" "data.frame"Que atributos tiene
attributes(Eje_cafetero_mun)$names
[1] "ISO" "NAME_0" "NAME_1" "NAME_2" "TYPE_2" "ENGTYPE_2"
[7] "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] 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
$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 identityst_crs(Eje_cafetero_mun)Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"Los límites
st_bbox(Eje_cafetero_mun) xmin ymin xmax ymax
-76.359497 4.135801 -74.648903 5.762600 Instalar paquetes:
#install.packages(c("lwgeom","units"))
library(lwgeom)
library(units)
library(mapview)###Calculo del área
(Eje_cafetero_mun$area <- st_area(Eje_cafetero_mun))Units: [m^2]
[1] 501640633 196738573 140637664 114270854 125392699 202221815 583686340
[8] 90164281 443872405 171875225 40294536 81859591 330148056 407560271
[15] 288383258 97649245 572255971 353539160 216927090 375746747 1033607646
[22] 111606211 528585935 440026533 105740896 113002677 35390260 207615048
[29] 62267516 84325761 143356378 215845105 84854008 119045635 188845404
[36] 83327252 392691604 255153726 148757309 175301329 70882318 101066578
[43] 73215138 30592240 164160228 566091856 651786111 929547212 149203658
[50] 510134104 212409204Para cambiar las unidades del área
set_units(Eje_cafetero_mun$area, km^2)Units: [km^2]
[1] 501.64063 196.73857 140.63766 114.27085 125.39270 202.22182 583.68634
[8] 90.16428 443.87241 171.87523 40.29454 81.85959 330.14806 407.56027
[15] 288.38326 97.64925 572.25597 353.53916 216.92709 375.74675 1033.60765
[22] 111.60621 528.58593 440.02653 105.74090 113.00268 35.39026 207.61505
[29] 62.26752 84.32576 143.35638 215.84510 84.85401 119.04564 188.84540
[36] 83.32725 392.69160 255.15373 148.75731 175.30133 70.88232 101.06658
[43] 73.21514 30.59224 164.16023 566.09186 651.78611 929.54721 149.20366
[50] 510.13410 212.40920library(readxl)
(my_data<- read_excel("CensoAgropecuario.xls"))###Union de datos del censo agropecuario al objeto sf
(Eje_cafetero_data<- my_data %>% filter(dep == c("Quindío", "Risaralda", "Caldas")))longer object length is not a multiple of shorter object lengthRenombrar los datos
(Eje_cafetero_data <- rename(Eje_cafetero_data, NAME_2 = mun))(Eje_cafetero_censo <- left_join(Eje_cafetero_mun, Eje_cafetero_data ))Joining, by = "NAME_2"Simple feature collection with 51 features and 14 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defsVisualización
suppressPackageStartupMessages(library(mapview))
Eje_cafetero_censo %>% mapview(zcol = "NAME_2", legend=TRUE, col.regions=sf.colors)##Ordenar los datos
Vamos a ordenar las municipales, basados en las Ha de bosque (primero mayor área de bosque):
(Eje_cafetero_clean <- Eje_cafetero_censo %>% arrange(desc(Ha_bosque_total)) %>% select (NAME_2, codigo_mun, Ha_bosque_total, Ha_agro_total, Ha_noagro_total, Ha_otros))Simple feature collection with 51 features and 6 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defsEje_cafetero_clean%>% mapview(zcol = "NAME_2", legend=TRUE, col.regions=sf.colors)