Puntos-Lineas-Poligonos en R
Luis Castillo
2022-08-22
INTRODUCCION
# To clean environment
rm(list = ls(all.names = TRUE))
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 507097 27.1 1126748 60.2 644200 34.5
## Vcells 900281 6.9 8388608 64.0 1634597 12.5library("sp")
library("rgdal")## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/luisc/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/luisc/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library("raster")
library("maptools")## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.library("tidyverse")## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()library ('gdalUtilities')
library ("landscapemetrics") library("sp")
library("rgdal")library("rgeos")## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-0
## Polygon checking: TRUElibrary("sf")## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE##
## Attaching package: 'sf'## The following object is masked from 'package:gdalUtilities':
##
## gdal_rasterizegetClass("Spatial") ## Class "Spatial" [package "sp"]
##
## Slots:
##
## Name: bbox proj4string
## Class: matrix CRS
##
## Known Subclasses:
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2Creación de objetos espaciales desde cero
Puntos espaciales
Vamos a empezar a generar un conjunto de 10 puntos en el cuadrado unitario [0, 1] y [0, 1]. La función runif() crea números aleatorios que siguen una distribución uniforme mientras que la función round() devuelve valores redondeados al número especificado de lugares decimales.
xc = round(runif(10), 2)
yc = round(runif(10), 2)
xy = cbind(xc, yc)
xy## xc yc
## [1,] 0.78 0.53
## [2,] 0.09 0.63
## [3,] 0.61 0.52
## [4,] 0.08 0.63
## [5,] 0.24 0.35
## [6,] 0.05 0.90
## [7,] 0.20 0.97
## [8,] 0.58 0.35
## [9,] 0.12 0.17
## [10,] 0.82 0.02class(xy)## [1] "matrix" "array"La matriz xy se puede convertir en un objeto SpatialPoints mediante la función SpatialPoints. Puede observar el cambio del tipo de objeto utilizando la función class().
xy.sp = SpatialPoints(xy)
class(xy.sp)## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"xy.sp## class : SpatialPoints
## features : 10
## extent : 0.05, 0.82, 0.02, 0.97 (xmin, xmax, ymin, ymax)
## crs : NAplot(xy.sp, pch = 2)
xy.cc = coordinates(xy.sp)
class(xy.cc)## [1] "matrix" "array"Puede obtener las coordenadas y el límite espacial de xy.sp con las funciones coordenadas() y bbox(), respectivamente. Este tipo de información puede ser útil cuando trabaja con polígonos creados por otras personas.
coordinates(xy.sp)## xc yc
## [1,] 0.78 0.53
## [2,] 0.09 0.63
## [3,] 0.61 0.52
## [4,] 0.08 0.63
## [5,] 0.24 0.35
## [6,] 0.05 0.90
## [7,] 0.20 0.97
## [8,] 0.58 0.35
## [9,] 0.12 0.17
## [10,] 0.82 0.02xy.cc = coordinates(xy.sp)#To create a data.frame with the coordinates (it is used a lot in spatial modeling)
class(xy.cc)## [1] "matrix" "array"bbox(xy.sp)## min max
## xc 0.05 0.82
## yc 0.02 0.97#To build the data.frame (df) with the attributes
df = data.frame(z1 = round(5 + rnorm(10), 2), z2 = 20:29)
df ## z1 z2
## 1 5.95 20
## 2 5.41 21
## 3 5.36 22
## 4 4.50 23
## 5 3.95 24
## 6 5.80 25
## 7 4.40 26
## 8 3.56 27
## 9 6.30 28
## 10 5.20 29#To join the SpatialPoint (x.ysp) to the data.frame (df)
xy.spdf = SpatialPointsDataFrame(xy.sp, df)
xy.spdf## class : SpatialPointsDataFrame
## features : 10
## extent : 0.05, 0.82, 0.02, 0.97 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 2
## names : z1, z2
## min values : 3.56, 20
## max values : 6.3, 29class(xy.spdf)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"Líneas espaciales
Las líneas espaciales se recuperan de fuentes externas con frecuencia; sin embargo, el siguiente ejemplo muestra cómo construir un objeto de clase SpatialLines desde cero. Primero creamos los puntos que seguirán las líneas. En segundo lugar, creamos las líneas usando estos puntos y la función Line(). En tercer lugar, ponemos estas líneas en una lista por separado; estas líneas no son objeto espacial en este punto. Finalmente fusionamos ambas líneas y damos propiedades espaciales con la función SpatialLines(). Tenga en cuenta que los objetos de líneas deben recibir valores de identificación de caracteres y que estos valores deben ser únicos para los objetos de líneas combinados en objetos SpatialLines.
# To create the points that the lines follow.
l1 = cbind(c(1,2,3),c(3,2,2))
l1a = cbind(l1[,1]+.05,l1[,2]+.05)
l2 = cbind(c(1,2,3),c(1,1.5,1))
# To create the lines.
Sl1 = Line(l1)
Sl1a = Line(l1a)
Sl2 = Line(l2)
# To put these lines in a list separately.
S1 = Lines(list(Sl1, Sl1a), ID="a")
S2 = Lines(list(Sl2), ID="b")
#
Sl = SpatialLines(list(S1,S2))
summary(Sl)## Object of class SpatialLines
## Coordinates:
## min max
## x 1 3.05
## y 1 3.05
## Is projected: NA
## proj4string : [NA]plot(Sl)
Podemos crear un objeto SpatialLines con atributos adjuntando una tabla de atributos (data.frame). En los siguientes códigos, la función sapply() aplica la función slot() a cada objeto de línea para generar una ID.
df_2 = data.frame(z = c(1,2), row.names=sapply(slot(Sl, "lines"), function(x) slot(x, "ID")))
Sldf = SpatialLinesDataFrame(Sl, data = df_2)
summary(Sldf)## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x 1 3.05
## y 1 3.05
## Is projected: NA
## proj4string : [NA]
## Data attributes:
## z
## Min. :1.00
## 1st Qu.:1.25
## Median :1.50
## Mean :1.50
## 3rd Qu.:1.75
## Max. :2.00#To verify ID (a and b) and the structure of Sldf.
str(Sldf)## Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
## ..@ data :'data.frame': 2 obs. of 1 variable:
## .. ..$ z: num [1:2] 1 2
## ..@ lines :List of 2
## .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
## .. .. .. ..@ Lines:List of 2
## .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
## .. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 1 2 3 3 2 2
## .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
## .. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 1.05 2.05 3.05 3.05 2.05 2.05
## .. .. .. ..@ ID : chr "a"
## .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
## .. .. .. ..@ Lines:List of 1
## .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
## .. .. .. .. .. .. ..@ coords: num [1:3, 1:2] 1 2 3 1 1.5 1
## .. .. .. ..@ ID : chr "b"
## ..@ bbox : num [1:2, 1:2] 1 1 3.05 3.05
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NA#To create no-spatial polygons
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
Sr1@coords #To see coordinates## [,1] [,2]
## [1,] 2 2
## [2,] 4 3
## [3,] 4 5
## [4,] 1 4
## [5,] 2 2#To put the polygons in a list and add IDs ("s1","s2","s3").
Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3, Sr4), "s3/s4")Tenga en cuenta que los objetos Polygons deben recibir valores de ID de carácter y que estos valores deben ser únicos para los objetos Polygons combinados en un objeto SpatialPolygons. Ahora vamos a crear el polígono espacial.
#To create the final spatial polygons and plot it.
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
plot(SpP, col = 1:3)
Como hicimos con las líneas y los puntos espaciales, podemos adjuntar atributos a los polígonos espaciales. Los polígonos con atributos (objetos de la clase SpatialPolygonsDataFrame) se construyen a partir del objeto SpatialPolygons (topología) y los atributos (data.frame). Los nombres de fila del marco de datos de atributos se comparan con las ranuras de ID del objeto SpatialPolygons, y las filas del marco de datos se reordenarán si es necesario.
attr = data.frame(a=1:3, row.names=c("s3/s4", "s2", "s1"))
attr## a
## s3/s4 1
## s2 2
## s1 3SrDf = SpatialPolygonsDataFrame(SpP, attr)
SrDf## class : SpatialPolygonsDataFrame
## features : 3
## extent : 1, 10, 2, 5 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 1
## names : a
## min values : 1
## max values : 3as(SrDf, "data.frame")## a
## s1 3
## s2 2
## s3/s4 1spplot(SrDf)# the function spplot() is a method for plotting spatial data with attributes.
Cargar y operar archivos de formas existentes en R
Los datos necesarios para visualizar un solo archivo de forma en Arc o QGIS se dividen en muchos archivos (p. ej., cpg, .dbf, .prj, etc.), el más conocido de los cuales es el archivo sagrado .SHP (¿conoce la función de los otros archivos?). Obtener toda esta información en un objeto espacial organizado en R es fácil y, de hecho, como suele ser el caso en R, existen varios paquetes y funciones que lo harán por usted. Estas son algunas de las opciones más utilizadas para leer .SHP de los paquetes rgdal, maptools, raster y sf:
library(rgdal)
DEPARTMENTS_1 = readOGR(dsn = "C:/Users/luisc/Downloads", layer = "MGN_DPTO_POLITICO") #no .shp required## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\luisc\Downloads", layer: "MGN_DPTO_POLITICO"
## with 33 features
## It has 10 fields
## Integer64 fields read as strings: DPTO_ANO_C DPTO_VGNClibrary(maptools)
DEPARTMENTS_2 = readShapePoly("C:/Users/luisc/Downloads/MGN_DPTO_POLITICO.shp")## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_readlibrary(raster)
DEPARTMENTS_3 = shapefile("C:/Users/luisc/Downloads/MGN_DPTO_POLITICO.shp")
library(sf)
DEPARTMENTS_4 = st_read("C:/Users/luisc/Downloads/MGN_DPTO_POLITICO.shp")## Reading layer `MGN_DPTO_POLITICO' from data source
## `C:\Users\luisc\Downloads\MGN_DPTO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## Geodetic CRS: WGS 84summary(DEPARTMENTS_1)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -81.735621 -66.84722
## y -4.229406 13.39473
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
## DPTO_CCDGO DPTO_CNMBR DPTO_ANO_C DPTO_ACT_A
## Length:33 Length:33 Length:33 Length:33
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## DPTO_NAREA DPTO_CSMBL DPTO_VGNC SHAPE_AREA
## Min. : 49.42 Length:33 Length:33 Min. :0.00405
## 1st Qu.: 20619.01 Class :character Class :character 1st Qu.:1.67949
## Median : 24139.40 Mode :character Mode :character Median :1.96503
## Mean : 34574.78 Mean :2.81313
## 3rd Qu.: 48353.22 3rd Qu.:3.93974
## Max. :109497.05 Max. :8.87748
##
## SHAPE_LEN AREA
## Min. : 0.6853 Min. : NA
## 1st Qu.: 9.5497 1st Qu.: NA
## Median :12.5785 Median : NA
## Mean :12.7108 Mean :NaN
## 3rd Qu.:17.2926 3rd Qu.: NA
## Max. :25.3560 Max. : NA
## NA's :33plot(DEPARTMENTS_1)
class(DEPARTMENTS_1@data)## [1] "data.frame"DEPARTMENTS_1@data## DPTO_CCDGO DPTO_CNMBR
## 0 05 ANTIOQUIA
## 1 23 CÓRDOBA
## 2 27 CHOCÓ
## 3 70 SUCRE
## 4 08 ATLÁNTICO
## 5 13 BOLÍVAR
## 6 47 MAGDALENA
## 7 20 CESAR
## 8 44 LA GUAJIRA
## 9 19 CAUCA
## 10 76 VALLE DEL CAUCA
## 11 41 HUILA
## 12 18 CAQUETÁ
## 13 50 META
## 14 15 BOYACÁ
## 15 25 CUNDINAMARCA
## 16 17 CALDAS
## 17 63 QUINDIO
## 18 66 RISARALDA
## 19 73 TOLIMA
## 20 52 NARIÑO
## 21 54 NORTE DE SANTANDER
## 22 68 SANTANDER
## 23 85 CASANARE
## 24 97 VAUPÉS
## 25 86 PUTUMAYO
## 26 94 GUAINÍA
## 27 99 VICHADA
## 28 91 AMAZONAS
## 29 95 GUAVIARE
## 30 88 ARCHIPIÉLAGO DE SAN ANDRÉS, PROVIDENCIA Y SANTA CATALINA
## 31 81 ARAUCA
## 32 11 BOGOTÁ, D.C.
## DPTO_ANO_C DPTO_ACT_A
## 0 1886 Constitucion Politica de 1886
## 1 1951 Ley 9 del 18 de Diciembre de 1951
## 2 1947 Ley 13 del 3 de Noviembre de 1947
## 3 1966 Ley 47 del 8 de Agosto de 1966
## 4 1910 Ley 21 de 1910
## 5 1886 Constitucion Politica de 1886
## 6 1964 1964
## 7 1967 Ley 25 21 de junio de 1967
## 8 1964 Acto Legislativo No. 1 de Diciembre 28 de 1964
## 9 1857 15 de junio de 1857
## 10 1910 Decreto No 340 de 16 de Abril de 1910
## 11 1905 Ley 46 de 1905
## 12 1981 Ley 78 del 29 de Diciembre de 1981
## 13 1959 Ley 118 del 16 de Diciembre de 1959
## 14 1886 Constitucion Politica de 1886
## 15 1886 Constitucion Politica de 1886
## 16 1905 11 de Abril de 1905
## 17 1966 Ley 2 TM de 1966
## 18 1966 Ley 70 del 1 de Diciembre de 1966
## 19 1909 Ley 65 de Noviembre de 1909
## 20 1904 Ley 1 de 1904
## 21 1910 Ley 25 de 1910
## 22 1910 Ley 25 14 de Julio de 1910
## 23 1991 5 de Julio Constitucion Politica de 1991
## 24 1991 Articulo 309 Constitucion Politica de 1991
## 25 1991 Articulo 309 Constitucion Politica de 1991
## 26 1991 Articulo 309 Constitucion Politica de 1991
## 27 1991 5 de Julio Constitucion Politica de 1991
## 28 1991 Dcto. 2274 del 4 de Octubre de la Constitución Política 1991
## 29 1991 5 de Julio Constitucion Politica de 1991
## 30 1991 Artículo 310 Constitucion Politica de 1991
## 31 1991 5 de Julio Constitucion Politica de 1991
## 32 1538 Constitucion Politica de 1886
## DPTO_NAREA DPTO_CSMBL DPTO_VGNC SHAPE_AREA SHAPE_LEN AREA
## 0 62804.71025 3 2020 5.134915412 21.4437936 NA
## 1 25086.54697 3 2020 2.057533219 9.6915160 NA
## 2 48353.21903 3 2020 3.939743093 20.6344626 NA
## 3 10591.84594 3 2020 0.870810349 8.5708693 NA
## 4 3313.81015 3 2020 0.273769829 2.5446507 NA
## 5 26719.21141 3 2020 2.195576862 16.2348170 NA
## 6 23135.93870 3 2020 1.909266273 10.8159579 NA
## 7 22565.30721 3 2020 1.858204407 12.5784592 NA
## 8 20619.00959 3 2020 1.706874490 10.7869049 NA
## 9 31242.91479 3 2020 2.534419222 13.9502629 NA
## 10 20665.54452 3 2020 1.679486591 12.6508699 NA
## 11 18141.66055 3 2020 1.474219177 10.3355646 NA
## 12 92831.12119 3 2020 7.540227855 21.2158385 NA
## 13 82799.17670 3 2020 6.733634872 18.1909537 NA
## 14 23138.04813 3 2020 1.888390832 15.9064914 NA
## 15 22370.48873 3 2020 1.823631133 13.1189510 NA
## 16 7425.22167 3 2020 0.605497805 6.6558442 NA
## 17 1933.57085 3 2020 0.157423729 2.5548254 NA
## 18 3556.77445 3 2020 0.289790632 4.8436149 NA
## 19 24139.40123 3 2020 1.965026841 9.5497261 NA
## 20 31497.57257 3 2020 2.548474605 12.8320974 NA
## 21 21856.75425 3 2020 1.792325441 10.7079284 NA
## 22 30561.51495 3 2020 2.499105532 11.8777470 NA
## 23 44394.23977 3 2020 3.615063148 12.1327536 NA
## 24 53299.28001 3 2020 4.313810177 20.1298337 NA
## 25 25976.28311 3 2020 2.107964833 12.7079225 NA
## 26 71289.35468 3 2020 5.747937395 21.1790507 NA
## 27 100063.37060 3 2020 8.100680492 17.2926130 NA
## 28 109497.05379 3 2020 8.877479989 25.3559774 NA
## 29 55575.23316 3 2020 4.511457243 19.3967890 NA
## 30 49.42425 3 2020 0.004049898 0.6852594 NA
## 31 23851.25706 3 2020 1.944156904 9.1245856 NA
## 32 1622.85260 3 2020 0.132207854 3.7604533 NADEPARTMENTS_1@proj4string## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## 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]]Algunas operaciones poligonales
Subconjunto
La primera operación poligonal principal es el subconjunto de atributos. Vamos a usar la función subset(). Recuerde agregar doble igual a "==" para seleccionar el departamento (META) de lo contrario la selección falla.
names(DEPARTMENTS_1)# To see the names of the attributes.## [1] "DPTO_CCDGO" "DPTO_CNMBR" "DPTO_ANO_C" "DPTO_ACT_A" "DPTO_NAREA"
## [6] "DPTO_CSMBL" "DPTO_VGNC" "SHAPE_AREA" "SHAPE_LEN" "AREA"head(DEPARTMENTS_1)# To see the first five rows## DPTO_CCDGO DPTO_CNMBR DPTO_ANO_C DPTO_ACT_A DPTO_NAREA
## 0 05 ANTIOQUIA 1886 Constitucion Politica de 1886 62804.71
## 1 23 CÓRDOBA 1951 Ley 9 del 18 de Diciembre de 1951 25086.55
## 2 27 CHOCÓ 1947 Ley 13 del 3 de Noviembre de 1947 48353.22
## 3 70 SUCRE 1966 Ley 47 del 8 de Agosto de 1966 10591.85
## 4 08 ATLÁNTICO 1910 Ley 21 de 1910 3313.81
## 5 13 BOLÍVAR 1886 Constitucion Politica de 1886 26719.21
## DPTO_CSMBL DPTO_VGNC SHAPE_AREA SHAPE_LEN AREA
## 0 3 2020 5.1349154 21.443794 NA
## 1 3 2020 2.0575332 9.691516 NA
## 2 3 2020 3.9397431 20.634463 NA
## 3 3 2020 0.8708103 8.570869 NA
## 4 3 2020 0.2737698 2.544651 NA
## 5 3 2020 2.1955769 16.234817 NAMETA = subset(DEPARTMENTS_1, DPTO_CNMBR == "META")
plot(META)
META## class : SpatialPolygonsDataFrame
## features : 1
## extent : -74.93335, -71.07753, 1.62084, 4.899101 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : DPTO_CCDGO, DPTO_CNMBR, DPTO_ANO_C, DPTO_ACT_A, DPTO_NAREA, DPTO_CSMBL, DPTO_VGNC, SHAPE_AREA, SHAPE_LEN, AREA
## value : 50, META, 1959, Ley 118 del 16 de Diciembre de 1959, 82799.176704, 3, 2020, 6.73363487204, 18.1909537045, NAMETA_HUILA = subset(DEPARTMENTS_1, DPTO_CNMBR == "META" | DPTO_CNMBR == "HUILA")
plot(META_HUILA)
META_HUILA## class : SpatialPolygonsDataFrame
## features : 2
## extent : -76.62466, -71.07753, 1.552125, 4.899101 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : DPTO_CCDGO, DPTO_CNMBR, DPTO_ANO_C, DPTO_ACT_A, DPTO_NAREA, DPTO_CSMBL, DPTO_VGNC, SHAPE_AREA, SHAPE_LEN, AREA
## min values : 41, HUILA, 1905, Ley 118 del 16 de Diciembre de 1959, 18141.6605485, 3, 2020, 1.47421917699, 10.3355645749, NA
## max values : 50, META, 1959, Ley 46 de 1905, 82799.176704, 3, 2020, 6.73363487204, 18.1909537045, NA#Group 1
Size = mean(DEPARTMENTS_1$DPTO_NAREA)#To obtain the mean
Size## [1] 34574.78GROUP_1 = subset(DEPARTMENTS_1, DPTO_NAREA >= mean(DEPARTMENTS_1$DPTO_NAREA) & DPTO_ANO_C <= 1950)
GROUP_1## class : SpatialPolygonsDataFrame
## features : 2
## extent : -77.88378, -73.88128, 3.964883, 8.873974 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : DPTO_CCDGO, DPTO_CNMBR, DPTO_ANO_C, DPTO_ACT_A, DPTO_NAREA, DPTO_CSMBL, DPTO_VGNC, SHAPE_AREA, SHAPE_LEN, AREA
## min values : 05, ANTIOQUIA, 1886, Constitucion Politica de 1886, 48353.2190281, 3, 2020, 3.93974309301, 20.634462642, NA
## max values : 27, CHOCÓ, 1947, Ley 13 del 3 de Noviembre de 1947, 62804.7102519, 3, 2020, 5.1349154118, 21.4437936323, NAplot(GROUP_1)
#Group 2
Qs = quantile(DEPARTMENTS_1$DPTO_NAREA) #To obtain the quartiles.
Qs## 0% 25% 50% 75% 100%
## 49.42425 20619.00959 24139.40123 48353.21903 109497.05379Qs[2]## 25%
## 20619.01GROUP_2 = subset(DEPARTMENTS_1, DPTO_NAREA <= Qs[2] & DPTO_ANO_C >= 1950)
GROUP_2## class : SpatialPolygonsDataFrame
## features : 5
## extent : -81.73562, -71.11296, 4.075311, 13.39473 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : DPTO_CCDGO, DPTO_CNMBR, DPTO_ANO_C, DPTO_ACT_A, DPTO_NAREA, DPTO_CSMBL, DPTO_VGNC, SHAPE_AREA, SHAPE_LEN, AREA
## min values : 44, ARCHIPIÉLAGO DE SAN ANDRÉS, PROVIDENCIA Y SANTA CATALINA, 1964, Acto Legislativo No. 1 de Diciembre 28 de 1964, 49.42425119, 3, 2020, 0.00404989764595, 0.685259427799, NA
## max values : 88, SUCRE, 1991, Ley 70 del 1 de Diciembre de 1966, 20619.0095911, 3, 2020, 1.70687449047, 10.7869048866, NAplot(GROUP_2)
Recortar (recortar)
También puede recortar funciones de entrada que se superponen a una función de clip. Puede usar recortes para cortar una parte de una clase de entidad usando una o más de las entidades en otra clase de entidad como cortador de galletas con la función crop(). Vamos a cargar los .SHP de la amazonia colombiana para recortar los departamentos que estan en la region.
library(rgeos)
AMAZONIA = readOGR(dsn = "C:/Users/luisc/Downloads", layer = "Lim_Amazonia")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\luisc\Downloads", layer: "Lim_Amazonia"
## with 1 features
## It has 5 fields
## Integer64 fields read as strings: objectidDEP_AMAZONIA = crop(x = DEPARTMENTS_1, y = AMAZONIA)
DEP_AMAZONIA## class : SpatialPolygonsDataFrame
## features : 11
## extent : -77.67061, -66.84722, -4.228294, 4.948186 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : DPTO_CCDGO, DPTO_CNMBR, DPTO_ANO_C, DPTO_ACT_A, DPTO_NAREA, DPTO_CSMBL, DPTO_VGNC, SHAPE_AREA, SHAPE_LEN, AREA
## min values : 18, AMAZONAS, 1857, 15 de junio de 1857, 18141.6605485, 3, 2020, 1.47421917699, 10.3355645749, NA
## max values : 99, VICHADA, 1991, Ley 78 del 29 de Diciembre de 1981, 109497.053788, 3, 2020, 8.87747998889, 25.3559773719, NAplot(DEP_AMAZONIA)
Eliminar (borrar)
Además, puedes eliminar polígonos aplicando la función erase(). Vamos a delatar el polígono META_HUILA (crear previamente).
COL_ERASE = erase(x = DEPARTMENTS_1, y = META_HUILA)
COL_ERASE## class : SpatialPolygonsDataFrame
## features : 31
## extent : -81.73562, -66.84722, -4.229406, 13.39473 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : DPTO_CCDGO, DPTO_CNMBR, DPTO_ANO_C, DPTO_ACT_A, DPTO_NAREA, DPTO_CSMBL, DPTO_VGNC, SHAPE_AREA, SHAPE_LEN, AREA
## min values : 05, AMAZONAS, 1538, 11 de Abril de 1905, 49.42425119, 3, 2020, 0.00404989764595, 0.685259427799, NA
## max values : 99, VICHADA, 1991, Ley 9 del 18 de Diciembre de 1951, 109497.053788, 3, 2020, 8.87747998889, 25.3559773719, NAplot(COL_ERASE, col ="blue")
Agregado (disolver)
Puede agregar polígonos en un objeto SpatialPolygons, de acuerdo con el vector de ID que especifica qué polígonos de entrada pertenecen a qué polígonos de salida; los límites internos se disuelven. Este tutorial presenta una opción para agregar polígonos usando la función gUnaryUnion() de la biblioteca rgeos. Necesitamos leer los municipios de Colombia (MUNI) en este punto del tutorial y crear un subconjunto de los municipios del departamento del Cauca.
MUNI = readOGR(dsn = "C:/Users/luisc/Downloads", layer = "MGN_MPIO_POLITICO")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\luisc\Downloads", layer: "MGN_MPIO_POLITICO"
## with 1121 features
## It has 12 fields
## Integer64 fields read as strings: MPIO_VGNCtable(MUNI$DPTO_CNMBR)##
## AMAZONAS
## 11
## ANTIOQUIA
## 125
## ARAUCA
## 7
## ARCHIPIÉLAGO DE SAN ANDRÉS, PROVIDENCIA Y SANTA CATALINA
## 2
## ATLÁNTICO
## 23
## BOGOTÁ, D.C.
## 1
## BOLÍVAR
## 46
## BOYACÁ
## 123
## CALDAS
## 27
## CAQUETÁ
## 16
## CASANARE
## 19
## CAUCA
## 42
## CESAR
## 25
## CHOCÓ
## 30
## CÓRDOBA
## 30
## CUNDINAMARCA
## 116
## GUAINÍA
## 8
## GUAVIARE
## 4
## HUILA
## 37
## LA GUAJIRA
## 15
## MAGDALENA
## 30
## META
## 29
## NARIÑO
## 64
## NORTE DE SANTANDER
## 40
## PUTUMAYO
## 13
## QUINDIO
## 12
## RISARALDA
## 14
## SANTANDER
## 87
## SUCRE
## 26
## TOLIMA
## 47
## VALLE DEL CAUCA
## 42
## VAUPÉS
## 6
## VICHADA
## 4CAUCA_MUNI = MUNI[MUNI$DPTO_CNMBR == "CAUCA",]
plot(CAUCA_MUNI,col = "green")
library(rgeos)
CAUCA_MUNI_dis_2 = gUnaryUnion(CAUCA_MUNI)
plot(CAUCA_MUNI_dis_2, col = "gold")
Proyecciones y centroides
La función gCentroid() calcula el centroide de la geometría dada. Vamos a calcularlos para CAUCA_MUNI
CAUCA_MUNI@proj4string## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## 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]]CENT = gCentroid(CAUCA_MUNI, byid=T)
plot(CENT) 
plot(CAUCA_MUNI)
plot(CENT, add=T, col = "red")
Ahora, podemos estimar la matriz de distancias entre los centroides. Sin embargo, R muestra un calentamiento: el objeto espacial 1 no se proyecta; GEOS espera coordenadas planas. Usted como analista espacial sabe que la estimación correcta de distancias o áreas necesita una proyección plana, y nuestros polígonos están en proyección geográfica.
dist_matrix <- gDistance(spgeom1 = CENT, spgeom2 = NULL, byid=T)## Warning in RGEOSDistanceFunc(spgeom1, spgeom2, byid, "rgeos_distance"): Spatial
## object 1 is not projected; GEOS expects planar coordinatesdist_matrix[1:4, 1:4] ## 299 300 301 302
## 299 0.0000000 0.1815861 1.530989 0.1026708
## 300 0.1815861 0.0000000 1.688305 0.1822681
## 301 1.5309886 1.6883047 0.000000 1.6083476
## 302 0.1026708 0.1822681 1.608348 0.0000000Por lo tanto, necesitamos cambiar la proyección geográfica de CAUCA_MUNI a una proyección plana. Vamos a usar UTM en la zona 18-norte, la UTM apropiada para el Cauca, usando la función spTransform. Consulte los argumentos del sistema de referencia de coordenadas o CRS, la nomenclatura de este argumento se explica aquí . ¿Cuáles son las unidades de las distancias?
CAUCA_MUNI_UTM18 = spTransform(CAUCA_MUNI, CRSobj="+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
proj4string(CAUCA_MUNI_UTM18)## Warning in proj4string(CAUCA_MUNI_UTM18): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html## [1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"CENT_UTM18 = gCentroid(CAUCA_MUNI_UTM18, byid=T)
plot(CENT_UTM18)
plot(CAUCA_MUNI_UTM18)
plot(CENT_UTM18, add=T, col = "purple")
dist_matrix_UTM18 <- gDistance(spgeom1 = CENT_UTM18, spgeom2 = NULL, byid=T)
dist_matrix_UTM18[1:4, 1:4]## 299 300 301 302
## 299 0.00 20155.21 169521.4 11355.53
## 300 20155.21 0.00 186985.0 20256.38
## 301 169521.40 186985.05 0.0 178051.07
## 302 11355.53 20256.38 178051.1 0.00AREA_UTM18 <- gArea(CAUCA_MUNI_UTM18, byid=T)
AREA_UTM18## 299 300 301 302 303 304 305
## 264286513 325348008 56525769 100473712 864837702 234025723 207923909
## 306 307 308 309 310 311 312
## 515015274 696381658 189353259 510430551 69557107 1790238326 753612832
## 313 314 315 316 317 318 319
## 1102765245 180917516 108920441 846767478 168745181 418184423 517210446
## 320 321 322 323 324 325 326
## 3612387039 679750940 514927923 391583600 135582723 200885003 491285730
## 327 328 329 330 331 332 520
## 403724395 2733595381 2561327119 3359786025 2052632936 81467763 479239409
## 521 522 523 524 525 526 527
## 236218575 773974389 412257773 796089724 434156019 551163152 354069693tampones
La zona de influencia es una técnica comúnmente utilizada para determinar áreas próximas a ciertos objetos. Tradicionalmente, los topes son circulares.
house_buffers_1 <- gBuffer(CENT_UTM18, width = 5000, byid = TRUE)
plot(house_buffers_1)
Pero R es tan flexible que también puede obtener buffers cuadrados y lineales. Sin embargo, el lineal es incompatible con las geometrías SpatialPoints.
house_buffers_3 <- gBuffer(CAUCA_MUNI_UTM18, width = 5000, byid = TRUE, capStyle = "FLAT")
plot(house_buffers_3)
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.