Review of Geometry and Topology
Ivan Lizarazo
August 23, 2016
Introduction
This is an example of testing functionalities of sp and rgeos libraries. The sp library is the basis for spatial classes in R. The rgeos library is the topology engine; it is based on the GEOS library which is a port of the Java Topology Suite (JTS). JTS and GEOS only handle planar geometry, so spatial objects should be projected. JTS and GEOS place coordinates on a very fine grid, scaling and rounding them, so that “very close” coordinates are placed on the same grid node. This leads to the concept of scaling factor which can be returned using the getScale() function.
JTS and GEOS follow the OpenGeospatial Consortium (OGC) Simple Features Specification. In this standard, polygons may have one and only one exterior boundary ring, and an unlimited number of interior boundaries (i.e. holes). This design may lead to several issues when working with polygonal objects with multiple exterior boundaries (i.e. a country made up of several islands) (Bivand et. al, 2013)
Simple Features Specification
Simple Features (officially Simple Feature Access) is both an OGC and ISO standard ISO 19125 that specifies a common storage and access model of mostly two-dimensional geographical data (point, line, polygon, multi-point, multi-line, etc.)
The ISO 19125 standard comes in two parts. Part one, ISO 19125-1 (SFA-CA for “common architecture”), defines a model for two-dimensional simple features, with linear interpolation between vertices. The data model defined in SFA-CA is a hierarchy of classes. This part also defines representation using Well-Known Text (and Binary). Part 2 of the standard, ISO 19125-2 (SFA-SQL), defines an implementation using SQL. The OpenGIS standard(s) cover implementations in CORBA and OLE/COM as well, although these have lagged behind the SQL one and are not standardized by ISO.
The ISO/IEC 13249-3 SQL/MM Spatial extends the Simple Features data model mainly with circular interpolations (e.g. circular arcs) and adds other features like coordinate transformations and methods for validating geometries as well as Geography Markup Language support.
Download data
Before starting, download following shapefiles provided by DANE and DIVA-GIS: (i) Marco geoestadistico Depto Cundinamarca (http://www.dane.gov.co/DescargaMGN/Departamento/MGN2005_25_CUNDINAMARCA.zip); (ii) centros poblados DANE 2012 (http://goo.gl/ehSSVo); (iii) Colombian inland water (DIVA-GIS); (iv) Colombian roads (DIVA-GIS); and (v) Colombian NBI (Necesidades Basicas Insatisfechas) DANE
Unzip the two shapefiles in your working directory.
Load objects of interest
As with any project, the starting point is to load the data we’ll be using and to understand its properties.
# define working directory
# setwd("./documentos/rpubs")
# load the packages we'll be using
library(rgdal)## Loading required package: sp## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3library(sp)
library(rgeos)## rgeos version: 0.3-19, (SVN revision 524)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-3
## Polygon checking: TRUE#
# list files:
list.files(path="./25_CUND/", pattern="shp")## [1] "MGN05_DPTO_25_CUN.shp" "MGN05_DPTO_25_CUN.shp.xml"
## [3] "MGN05_MPIO_25_CUN.shp" "MGN05_MPIO_25_CUN.shp.xml"
## [5] "MGN05_MZ_25CUN.shp" "MGN05_MZ_25CUN.shp.xml"
## [7] "MGN05_PERCEN_25_CUN.shp" "MGN05_PERCEN_25_CUN.shp.xml"
## [9] "MGN05_SCCR_25_CUN.shp" "MGN05_SCCR_25_CUN.shp.xml"
## [11] "MGN05_SCCU_25CUN.shp" "MGN05_SCCU_25CUN.shp.xml"
## [13] "MGN05_SCTR_25_CUN.shp" "MGN05_SCTR_25_CUN.shp.xml"
## [15] "MGN05_SCTU_25_CUN.shp" "MGN05_SCTU_25_CUN.shp.xml"# note misspelled names !!!
# select data:
# let's choose a shapefile (spatial data) and a excel file (non spatial)
#
# read municipios data:
mun <- readOGR('./25_CUND/MGN05_MPIO_25_CUN.shp', layer='MGN05_MPIO_25_CUN')## OGR data source with driver: ESRI Shapefile
## Source: "./25_CUND/MGN05_MPIO_25_CUN.shp", layer: "MGN05_MPIO_25_CUN"
## with 116 features
## It has 5 fields#
summary(mun)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -74.890282 -73.051163
## y 3.727698 5.837239
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CNMBR MPIO_CCNCT
## Min. :436.0 25:116 001 : 1 AGUA DE DIOS: 1 25001 : 1
## 1st Qu.:616.8 019 : 1 ALBAN : 1 25019 : 1
## Median :657.5 035 : 1 ANAPOIMA : 1 25035 : 1
## Mean :665.1 040 : 1 ANOLAIMA : 1 25040 : 1
## 3rd Qu.:714.0 053 : 1 APULO : 1 25053 : 1
## Max. :790.0 086 : 1 ARBELAEZ : 1 25086 : 1
## (Other):110 (Other) :110 (Other):110#
# what is the CRS of these objects
mun$proj4string## NULL# plot municipios
# start by plotting municipios
plot(mun, border="wheat3", col="wheat1", axes=TRUE)
# add the river lines
# lines(rios, col="dodgerblue3")
# add municipios labels (using trial and error for placement)
text(coordinates(mun), labels = mun$MPIO_CNMBR, col="darkred",
cex=0.6, font=2, offset=0.5, adj=c(0,2))
# add a plot title and legend
title("Municipios de Cundinamarca")
#
# read NBI data:
csv.file <- "./25_CUND/nbi2.csv"
NBI <- read.csv(file=csv.file, header=TRUE, sep=",")
head(NBI)## CODIGO_DEPTO NOMBRE_DEPTO CODIGO_MUN NOMBRE_MUN NBI_CABEC NBI_RESTO
## 1 25 CUNDINAMARCA 1 AGUA DE DIOS 15.17 38.92
## 2 25 CUNDINAMARCA 19 ALBAN 20.57 27.62
## 3 25 CUNDINAMARCA 35 ANAPOIMA 21.25 36.42
## 4 25 CUNDINAMARCA 40 ANOLAIMA 16.65 31.20
## 5 25 CUNDINAMARCA 53 ARBELAEZ 20.66 28.30
## 6 25 CUNDINAMARCA 86 BELTRAN 33.23 33.99
## NBI_TOTAL
## 1 20.672
## 2 25.783
## 3 29.981
## 4 26.859
## 5 25.246
## 6 33.860summary(NBI)## CODIGO_DEPTO NOMBRE_DEPTO CODIGO_MUN NOMBRE_MUN
## Min. :25 CUNDINAMARCA:116 Min. : 1.0 AGUA DE DIOS: 1
## 1st Qu.:25 1st Qu.:287.5 ALBAN : 1
## Median :25 Median :478.0 ANAPOIMA : 1
## Mean :25 Mean :485.9 ANOLAIMA : 1
## 3rd Qu.:25 3rd Qu.:755.0 APULO : 1
## Max. :25 Max. :899.0 ARBELAEZ : 1
## (Other) :110
## NBI_CABEC NBI_RESTO NBI_TOTAL
## Min. : 4.45 Min. : 9.35 Min. : 7.112
## 1st Qu.:13.46 1st Qu.:26.15 1st Qu.:21.408
## Median :18.22 Median :34.97 Median :29.549
## Mean :19.43 Mean :35.25 Mean :30.230
## 3rd Qu.:24.50 3rd Qu.:43.64 3rd Qu.:37.899
## Max. :41.88 Max. :76.40 Max. :68.487
## Join attributes
Joining attributes to spatial objects is essential for any analysis. Here, we’ll join “unsatisfied basic needs” (NBI) to Cundinamarca municipalities.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionmun$NOMBRE_MUN <- mun$MPIO_CNMBR
mun@data <- left_join(mun@data, NBI)## Joining, by = "NOMBRE_MUN"## Warning in left_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining
## factors with different levels, coercing to character vectorhead(mun@data) # the new data has been added## OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CNMBR MPIO_CCNCT NOMBRE_MUN
## 1 436 25 372 JUNIN 25372 JUNIN
## 2 437 25 317 GUACHETA 25317 GUACHETA
## 3 438 25 377 LA CALERA 25377 LA CALERA
## 4 439 25 386 LA MESA 25386 LA MESA
## 5 440 25 394 LA PALMA 25394 LA PALMA
## 6 592 25 398 LA PEÑA 25398 LA PEÑA
## CODIGO_DEPTO NOMBRE_DEPTO CODIGO_MUN NBI_CABEC NBI_RESTO NBI_TOTAL
## 1 25 CUNDINAMARCA 372 9.29 36.77 34.108
## 2 25 CUNDINAMARCA 317 28.39 34.86 32.823
## 3 25 CUNDINAMARCA 377 4.45 19.33 13.358
## 4 25 CUNDINAMARCA 386 11.76 36.83 23.687
## 5 25 CUNDINAMARCA 394 16.64 60.41 42.753
## 6 NA <NA> NA NA NA NAmun$MPIO_CNMBR <- NULL
#mun$NBI_TOTAL <- as.double(mun$NBI_TOTAL)
#mun$NBI_CABEC <- as.double(mun$NBI_CABEC)
#mun$NBI_RESTO <- as.double(mun$NBI_RESTO)
summary(mun)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -74.890282 -73.051163
## y 3.727698 5.837239
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CCNCT NOMBRE_MUN
## Min. :436.0 25:116 001 : 1 25001 : 1 Length:116
## 1st Qu.:616.8 019 : 1 25019 : 1 Class :character
## Median :657.5 035 : 1 25035 : 1 Mode :character
## Mean :665.1 040 : 1 25040 : 1
## 3rd Qu.:714.0 053 : 1 25053 : 1
## Max. :790.0 086 : 1 25086 : 1
## (Other):110 (Other):110
## CODIGO_DEPTO NOMBRE_DEPTO CODIGO_MUN NBI_CABEC
## Min. :25 CUNDINAMARCA:111 Min. : 1.0 Min. : 4.45
## 1st Qu.:25 NA's : 5 1st Qu.:287.0 1st Qu.:13.41
## Median :25 Median :473.0 Median :18.03
## Mean :25 Mean :484.1 Mean :19.20
## 3rd Qu.:25 3rd Qu.:756.0 3rd Qu.:24.38
## Max. :25 Max. :899.0 Max. :37.18
## NA's :5 NA's :5 NA's :5
## NBI_RESTO NBI_TOTAL
## Min. : 9.35 Min. : 7.112
## 1st Qu.:25.82 1st Qu.:21.384
## Median :34.86 Median :29.256
## Mean :34.74 Mean :29.709
## 3rd Qu.:42.67 3rd Qu.:37.602
## Max. :72.95 Max. :64.726
## NA's :5 NA's :5Find municipios with high poverty
Let’s find municipios with a cabecera municipal having a NBI higher than 55%.
mun@data <- na.omit(mun@data)
high_nbi <- mun[mun$NBI_TOTAL > 55.0,]
head(high_nbi@data)## OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CCNCT NOMBRE_MUN CODIGO_DEPTO
## 21 608 25 518 25518 PAIME 25
## 101 775 25 823 25823 TOPAIPI 25
## 113 787 25 885 25885 YACOPI 25
## NOMBRE_DEPTO CODIGO_MUN NBI_CABEC NBI_RESTO NBI_TOTAL
## 21 CUNDINAMARCA 518 16.07 64.34 59.719
## 101 CUNDINAMARCA 823 29.49 69.61 63.581
## 113 CUNDINAMARCA 885 32.82 72.95 64.726plot(high_nbi, col="lightblue", axes=TRUE)
text(coordinates(high_nbi), labels = high_nbi$NOMBRE_MUN, col="darkgrey",
cex=1.0, font=2, offset=0.5, adj=c(0,2))
Read centros poblados
As you noticed when dowloading Centros poblados, DANE provides this dataset in different formats. It includes population centers for the whole Colombian territory. Let’s read and examine the corresponding shapefile.
# read centros poblados data
poblados <- readOGR('./25_CUND/poblaciones/poblaciones_dane_2012.shp', layer="poblaciones_dane_2012")## OGR data source with driver: ESRI Shapefile
## Source: "./25_CUND/poblaciones/poblaciones_dane_2012.shp", layer: "poblaciones_dane_2012"
## with 4843 features
## It has 11 fieldssummary(poblados)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -81.725786 -67.06488
## y -4.228267 13.38684
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## OBJECTID CPOB_CCDGO NOMBRE CODANEPOB
## Min. : 1 000 :1115 Pueblo Nuevo: 21 05001000: 2
## 1st Qu.:1212 001 : 448 San Isidro : 19 05147000: 2
## Median :2422 002 : 373 San Pedro : 18 05212000: 2
## Mean :2422 003 : 308 San Antonio : 17 05266000: 2
## 3rd Qu.:3632 004 : 248 San José : 16 05380000: 2
## Max. :4843 005 : 241 Santa Rosa : 15 13647002: 2
## (Other):2110 (Other) :4737 (Other) :4831
## CODANE Tipo Nom_DANE Shape_Leng
## 52835 : 131 C :1755 PUEBLO NUEVO: 21 Min. :0.001498
## 76109 : 79 CM :1104 SAN ISIDRO : 19 1st Qu.:0.010610
## 66001 : 67 CAS : 961 SAN JOSÉ : 18 Median :0.018074
## 23001 : 51 CP : 461 SAN ANTONIO : 17 Mean :0.033240
## 52001 : 41 IPD : 264 SAN PEDRO : 17 3rd Qu.:0.033025
## 76520 : 40 IP : 146 SANTA ROSA : 15 Max. :2.219656
## (Other):4434 (Other): 152 (Other) :4736
## Shape_Area LAT LON
## Min. :1.100e-07 Min. :-4.202 Min. :-81.72
## 1st Qu.:3.120e-06 1st Qu.: 3.847 1st Qu.:-76.13
## Median :8.510e-06 Median : 5.587 Median :-75.27
## Mean :7.173e-05 Mean : 6.008 Mean :-75.14
## 3rd Qu.:2.566e-05 3rd Qu.: 8.716 3rd Qu.:-74.12
## Max. :3.385e-02 Max. :13.384 Max. : 0.00
## # what is the CRS of these objects?
#
# plot municipios de Cundinamarca and a subset of poblados
plot(mun, col=c("lightgrey","lightblue"), axes=TRUE)
plot(poblados, add=TRUE)
text(coordinates(mun), labels = mun$NOMBRE_MUN, col="lightgrey",
cex=1.0, font=2, offset=0.5, adj=c(0,2))
Read roads and water
DIVA-GIS provides hydrographic datasets for any country. We’ll use here the water lines (i.e. streams) as well as water areas (i.e. rivers and lakes) for the whole Colombian territory.
roads <- readOGR('./COL_rds/COL_roads.shp', layer="COL_roads")## OGR data source with driver: ESRI Shapefile
## Source: "./COL_rds/COL_roads.shp", layer: "COL_roads"
## with 2593 features
## It has 5 fieldsrivers <- readOGR('./COL_wat/COL_water_areas_dcw.shp', layer="COL_water_areas_dcw")## OGR data source with driver: ESRI Shapefile
## Source: "./COL_wat/COL_water_areas_dcw.shp", layer: "COL_water_areas_dcw"
## with 1166 features
## It has 5 fieldsstreams <- readOGR('./COL_wat/COL_water_lines_dcw.shp', layer="COL_water_lines_dcw")## OGR data source with driver: ESRI Shapefile
## Source: "./COL_wat/COL_water_lines_dcw.shp", layer: "COL_water_lines_dcw"
## with 3673 features
## It has 5 fieldsplot(poblados[poblados$NOMBRE == "Barranquilla",], col=c("lightgrey","lightblue"), axes=TRUE)
plot(streams,col="darkblue", add=TRUE)
plot(rivers,col="darkblue", add=TRUE)
plot(roads,col="darkred", add=TRUE)
text(coordinates(poblados), labels = poblados$NOMBRE, col="darkgrey",
cex=2.0, font=2, offset=0.5, adj=c(0,2))
Project data
In order to use rgeos, a plane coordinate system is needed. For reprojecting data, the spTransform() function can be used. Once this reprojection is done, you can test spatial relationships betwen objects.
mun2 <- readOGR('./25_CUND/MGN05_MPIO_25_CUN.shp', layer='MGN05_MPIO_25_CUN')## OGR data source with driver: ESRI Shapefile
## Source: "./25_CUND/MGN05_MPIO_25_CUN.shp", layer: "MGN05_MPIO_25_CUN"
## with 116 features
## It has 5 fieldsmun_bog <- spTransform(mun2, CRS("+init=epsg:3116"))
poblados_bog <- spTransform(poblados, CRS("+init=epsg:3116"))
roads_bog <- spTransform(roads, CRS("+init=epsg:3116"))
rivers_bog <- spTransform(rivers, CRS("+init=epsg:3116"))
streams_bog <- spTransform(streams, CRS("+init=epsg:3116"))
summary(mun_bog)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 909771.1 1113871
## y 903969.8 1137245
## Is projected: TRUE
## proj4string :
## [+init=epsg:3116 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CNMBR MPIO_CCNCT
## Min. :436.0 25:116 001 : 1 AGUA DE DIOS: 1 25001 : 1
## 1st Qu.:616.8 019 : 1 ALBAN : 1 25019 : 1
## Median :657.5 035 : 1 ANAPOIMA : 1 25035 : 1
## Mean :665.1 040 : 1 ANOLAIMA : 1 25040 : 1
## 3rd Qu.:714.0 053 : 1 APULO : 1 25053 : 1
## Max. :790.0 086 : 1 ARBELAEZ : 1 25086 : 1
## (Other):110 (Other) :110 (Other):110summary(poblados_bog)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 166577.97 1782176
## y 23000.04 1983987
## Is projected: TRUE
## proj4string :
## [+init=epsg:3116 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## OBJECTID CPOB_CCDGO NOMBRE CODANEPOB
## Min. : 1 000 :1115 Pueblo Nuevo: 21 05001000: 2
## 1st Qu.:1212 001 : 448 San Isidro : 19 05147000: 2
## Median :2422 002 : 373 San Pedro : 18 05212000: 2
## Mean :2422 003 : 308 San Antonio : 17 05266000: 2
## 3rd Qu.:3632 004 : 248 San José : 16 05380000: 2
## Max. :4843 005 : 241 Santa Rosa : 15 13647002: 2
## (Other):2110 (Other) :4737 (Other) :4831
## CODANE Tipo Nom_DANE Shape_Leng
## 52835 : 131 C :1755 PUEBLO NUEVO: 21 Min. :0.001498
## 76109 : 79 CM :1104 SAN ISIDRO : 19 1st Qu.:0.010610
## 66001 : 67 CAS : 961 SAN JOSÉ : 18 Median :0.018074
## 23001 : 51 CP : 461 SAN ANTONIO : 17 Mean :0.033240
## 52001 : 41 IPD : 264 SAN PEDRO : 17 3rd Qu.:0.033025
## 76520 : 40 IP : 146 SANTA ROSA : 15 Max. :2.219656
## (Other):4434 (Other): 152 (Other) :4736
## Shape_Area LAT LON
## Min. :1.100e-07 Min. :-4.202 Min. :-81.72
## 1st Qu.:3.120e-06 1st Qu.: 3.847 1st Qu.:-76.13
## Median :8.510e-06 Median : 5.587 Median :-75.27
## Mean :7.173e-05 Mean : 6.008 Mean :-75.14
## 3rd Qu.:2.566e-05 3rd Qu.: 8.716 3rd Qu.:-74.12
## Max. :3.385e-02 Max. :13.384 Max. : 0.00
## summary(rivers_bog)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 448009.63 1803730
## y 23358.58 1869402
## Is projected: TRUE
## proj4string :
## [+init=epsg:3116 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## ISO COUNTRY F_CODE_DES
## COL:1166 Colombia:1166 Inland Water :1162
## Land Subject to Inundation: 4
##
##
##
##
##
## HYC_DESCRI NAME
## Non-Perennial/Intermittent/Fluctuating: 9 UNK :921
## Perennial/Permanent :1157 RIO CAQUETA : 10
## CIENAGA PONEDERA: 6
## RIO APAPORIS : 6
## RIO PUTUMAYO : 6
## (Other) :213
## NA's : 4summary(streams_bog)## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x 457691.00 1800970
## y 34364.58 1862950
## Is projected: TRUE
## proj4string :
## [+init=epsg:3116 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## F_CODE_DES HYC_DESCRI
## River/Stream:3673 Non-Perennial/Intermittent/Fluctuating: 86
## Perennial/Permanent :3587
##
##
##
##
##
## NAM ISO NAME_0
## UNK :3129 COL:3673 Colombia:3673
## RIO CAUCA : 23
## RIO CESAR : 20
## RIO SAN JORGE: 20
## RIO TOMO : 18
## RIO APAPORIS : 15
## (Other) : 448summary(roads_bog)## Object of class SpatialLinesDataFrame
## Coordinates:
## min max
## x 473784.1 1731413
## y 141999.9 1869917
## Is projected: TRUE
## proj4string :
## [+init=epsg:3116 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
## MED_DESCRI RTT_DESCRI F_CODE_DES ISO
## Unknown : 68 Primary Route : 306 Road :2589 COL:2593
## With Median : 1 Secondary Route:2215 Trail: 4
## Without Median:2520 Unknown : 68
## NA's : 4 NA's : 4
## ISOCOUNTRY
## COLOMBIA:2589
## NA's : 4
##
## # plot
plot(poblados_bog[poblados_bog$NOMBRE == "Puerto Wilches",], col=c("lightgrey","lightblue"), axes=TRUE)
plot(streams_bog, col="darkblue", add=TRUE)
plot(rivers_bog, col="darkblue", add=TRUE)
plot(roads_bog, col="darkred", add=TRUE)
text(coordinates(poblados_bog), labels = poblados_bog$NOMBRE, col="darkgrey",
cex=2.0, font=2, offset=0.5, adj=c(0,2))
Find roads intersecting Cundinamarca
gUnaryUnion() is the R equivalent to ArcGIS dissolve operation.
cundi <- gUnaryUnion(mun_bog, id = mun_bog@data$DPTO_DPTO_)
roads_cundi = gIntersection(roads_bog,cundi)
plot(cundi, border="lightgrey", axes=TRUE)
plot(poblados_bog, col="lightblue", add=TRUE)
plot(roads_cundi, col="darkred", lwd=2, add=TRUE)
Find streams intersecting Cundinamarca
gIntersection() is a quick way to intersect polygons using rgeos. Unfortunately, it doesn’t retrieve the associated dataframes and their attributes. Another option could be using the insersection() function available in raster package. A third way could be using the getIntersection() function in RFigisGeo.
streams_cundi = gIntersection(streams_bog,cundi)
plot(cundi, border="lightgrey", axes=TRUE)
plot(mun_bog, col="lightblue", add=TRUE)
plot(streams_cundi, col="darkblue", lwd=2, add=TRUE)
Find roads within 50 km of Girardot’s population center
rgeos provides a gBuffer() function appropriate for this task.
girardot <- poblados_bog[poblados_bog$NOMBRE == "Girardot",]
gir_buf = gBuffer(girardot, width=50000)
gir_buf## An object of class "SpatialPolygons"
## Slot "polygons":
## [[1]]
## An object of class "Polygons"
## Slot "Polygons":
## [[1]]
## An object of class "Polygon"
## Slot "labpt":
## [1] 920021.0 968961.2
##
## Slot "area":
## [1] 8914307557
##
## Slot "hole":
## [1] FALSE
##
## Slot "ringDir":
## [1] 1
##
## Slot "coords":
## x y
## [1,] 937292.9 919058.3
## [2,] 926116.2 916148.0
## [3,] 926093.3 916144.8
## [4,] 923022.8 915821.7
## [5,] 922991.6 915819.4
## [6,] 918424.0 915690.3
## [7,] 918388.3 915691.0
## [8,] 915721.6 915808.8
## [9,] 915658.5 915813.2
## [10,] 912651.1 916118.3
## [11,] 912584.8 916127.0
## [12,] 912521.8 916135.4
## [13,] 909992.5 916869.4
## [14,] 901027.2 918540.3
## [15,] 888016.7 925288.0
## [16,] 877521.7 935518.3
## [17,] 877514.7 935527.5
## [18,] 871710.1 945267.3
## [19,] 871239.5 946310.4
## [20,] 871521.8 946437.8
## [21,] 870271.9 948511.3
## [22,] 866826.2 964726.9
## [23,] 868047.9 974551.2
## [24,] 867743.0 974595.0
## [25,] 867821.2 975139.4
## [26,] 868361.3 977071.3
## [27,] 868871.9 981177.9
## [28,] 869125.4 981693.5
## [29,] 869801.2 986552.1
## [30,] 876632.0 1000614.3
## [31,] 876644.3 1000631.9
## [32,] 879206.0 1003991.7
## [33,] 879259.7 1004056.3
## [34,] 882413.7 1007507.4
## [35,] 882544.1 1007637.3
## [36,] 888486.6 1012698.1
## [37,] 888505.6 1012711.9
## [38,] 890980.1 1014395.2
## [39,] 891066.7 1014450.4
## [40,] 896126.6 1017277.1
## [41,] 896142.3 1017284.7
## [42,] 910061.0 1021666.8
## [43,] 910887.9 1021676.8
## [44,] 911199.7 1021777.8
## [45,] 920682.0 1022021.4
## [46,] 921275.0 1022082.3
## [47,] 921496.8 1022042.3
## [48,] 925523.4 1022145.7
## [49,] 934972.8 1019612.9
## [50,] 935043.0 1019600.3
## [51,] 935070.2 1019586.8
## [52,] 939363.3 1018436.2
## [53,] 939610.9 1018329.8
## [54,] 939900.5 1018204.3
## [55,] 939943.0 1018185.8
## [56,] 944792.7 1015029.0
## [57,] 945455.1 1014699.4
## [58,] 945460.5 1014696.6
## [59,] 946620.0 1014067.4
## [60,] 947234.6 1013723.3
## [61,] 947226.1 1013708.0
## [62,] 950776.8 1011866.9
## [63,] 950798.4 1011852.6
## [64,] 950985.9 1011727.8
## [65,] 951785.6 1011192.9
## [66,] 953988.1 1009632.1
## [67,] 954023.4 1009605.7
## [68,] 955910.4 1008118.2
## [69,] 955930.0 1008102.1
## [70,] 965485.3 997563.3
## [71,] 965490.9 997554.9
## [72,] 966540.0 994892.3
## [73,] 971758.8 984061.2
## [74,] 971765.8 984039.7
## [75,] 972440.3 979917.8
## [76,] 973434.2 977395.3
## [77,] 973476.2 977130.0
## [78,] 973933.3 973278.0
## [79,] 973979.4 972697.8
## [80,] 973626.3 972669.7
## [81,] 974190.3 969223.2
## [82,] 974190.5 969205.0
## [83,] 973236.1 958891.5
## [84,] 973230.2 958861.9
## [85,] 970868.2 950685.8
## [86,] 970853.9 950648.6
## [87,] 970852.0 950643.7
## [88,] 967730.6 944011.4
## [89,] 967711.1 943976.7
## [90,] 964347.0 940293.6
## [91,] 964006.9 939750.8
## [92,] 964146.0 939654.0
## [93,] 963702.2 939016.0
## [94,] 963593.8 939091.4
## [95,] 962324.8 937065.9
## [96,] 962302.0 937036.0
## [97,] 952862.3 927719.4
## [98,] 952540.1 927366.7
## [99,] 952521.9 927354.1
## [100,] 952464.2 927326.6
## [101,] 952334.0 927198.0
## [102,] 952310.8 927180.8
## [103,] 949124.4 925329.2
## [104,] 948428.1 924771.9
## [105,] 948403.5 924755.1
## [106,] 944612.0 922403.7
## [107,] 944547.5 922367.6
## [108,] 940820.7 920480.6
## [109,] 940700.0 920425.7
## [110,] 939748.8 920004.9
## [111,] 939707.3 919987.0
## [112,] 937813.1 919218.0
## [113,] 937306.0 919024.1
## [114,] 937292.9 919058.3
##
##
##
## Slot "plotOrder":
## [1] 1
##
## Slot "labpt":
## [1] 920021.0 968961.2
##
## Slot "ID":
## [1] "buffer"
##
## Slot "area":
## [1] 8914307557
##
##
##
## Slot "plotOrder":
## [1] 1
##
## Slot "bbox":
## min max
## x 866826.2 974190.5
## y 915690.3 1022145.7
##
## Slot "proj4string":
## CRS arguments:
## +init=epsg:3116 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defsout <- gIntersection(gir_buf, roads_bog)
plot(gir_buf, axes=TRUE)
plot(mun_bog, col="lightgray", add=TRUE)
plot(girardot, col="darkgreen", add=TRUE)
plot(out, col=adjustcolor("darkred",alpha=0.5), add=TRUE, border=rgb(0,0,0,alpha=0))
text(coordinates(girardot), labels = girardot$NOMBRE, col="darkgrey",
cex=2.0, font=2, offset=0.5, adj=c(0,2))
Now, it’s your turn to code
Do not forget to replicate this exercise. Furthermore, use your imagination (and the rgeos library help) to write your own code for testing various rgeos functions such as gArea, gBoundary, gBuffer, gCentroid, gContains, gConvexHull, gCrosses, gDifference, gDistance, gEnvelope, gEquals, gIntersects, gIsEmpty, gIsRing, gIsSimple, gIsValid, and gLength.
Remember: Do not copy & paste. If you are serious about learning R, type code by yourself.
REFERENCES
Bivand, R.S., Pebesma, E., Gomez-Rubio, V., 2013. Applied Spatial Data Analysis with R. Second Edition. Springer. 405 pp.
Bivand, R.S., Rundel, C., Pebesma, E., Hufthammer, K.O., 2016. Interface to Geometry Engine - Open Source (GEOS). Available at: https://cran.r-project.org/web/packages/rgeos/rgeos.pdf
Oyana T.J. & Margai, F.M., 2016. Spatial Analysis – Statistics, Visualization, and Computational Methods. CRC Press. 294 pp.