inst/materials/lectures/day2.md
In thk686/keitt.ssi.2019: Learning materials for Geospatial Data Analysis in R, part of the UT Summer Statistics Institute
Geospatial Data Analysis in R
author: Timothy H. Keitt
date: May 12, 2014
width: 1440
height: 900
Working with vector data
type: section
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Spatial layers
GIS building blocks
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
OGC Simple Feature Hierarchy
Working with vector data
OGC Simple Feature Hierarchy
Working with vector data
Complex versus simple features
Most spatial operators in R require simple features
Working with vector data
Complex versus simple features
Most spatial operators in R require simple features
Working with vector data
OGC Simple Features Well-Known-Text
Working with vector data
OGC Simple Features Well-Known-Text
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
RGDAL package
- OGR component of GDAL handles vector data (http://www.gdal.org/)
- Data source -> data layer -> simple features
- Written in C++
- Like GDAL, large collection of runtime drivers
- Encapsulated in
rgdal and experimental rgdal2 packages
- With
rgdal package, main use is to read data into sp classes
Working with vector data
RGDAL package
pacman::p_load(rgdal)
as.character(ogrDrivers()[[1]])
[1] "AeronavFAA" "AmigoCloud" "ARCGEN" "AVCBin"
[5] "AVCE00" "BNA" "Carto" "Cloudant"
[9] "CouchDB" "CSV" "CSW" "DGN"
[13] "DXF" "EDIGEO" "ElasticSearch" "ESRI Shapefile"
[17] "Geoconcept" "GeoJSON" "GeoRSS" "GFT"
[21] "GML" "GPKG" "GPSBabel" "GPSTrackMaker"
[25] "GPX" "HTF" "HTTP" "Idrisi"
[29] "JML" "KML" "MapInfo File" "Memory"
[33] "netCDF" "ODS" "OGR_GMT" "OGR_PDS"
[37] "OGR_SDTS" "OGR_VRT" "OpenAir" "OpenFileGDB"
[41] "OSM" "PCIDSK" "PDF" "PGDUMP"
[45] "PLSCENES" "REC" "S57" "SEGUKOOA"
[49] "SEGY" "Selafin" "SQLite" "SUA"
[53] "SVG" "SXF" "TIGER" "UK .NTF"
[57] "VDV" "VFK" "WAsP" "WFS"
[61] "XLSX" "XPlane"
Uniform access to many data sources
Uniform access to many data sources
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Vector data IO
- main functions
- readOGR
- writeOGR
vecdat = readOGR("example-data/continents", "continent")
OGR data source with driver: ESRI Shapefile
Source: "/Users/tkeitt/Dropbox/R/keitt.ssi.2019/inst/materials/lectures/example-data/continents", layer: "continent"
with 8 features
It has 1 fields
First argument (dsn) is the data source. The second argument is the layer within the data source.
Working with vector data
Vector data IO
An ESRI shapefile is often several files stored together in a directory. The directory is the data source and the .shp file is the layer. Other files hold spatial reference system and tabular field data.
dir("example-data") # see the data source
[1] "continents" "NCEAS sample" "ozone.gml" "ozone.xsd"
[5] "rainforest" "rainforest.tif"
dir("example-data/continents") # see the data layers
[1] "continent.dbf" "continent.prj" "continent.shp" "continent.shx"
Other formats are single file, in which case the data source is the file name and the layer is named within the file, typically the same as the file name without the extension.
Working with vector data
Vector data IO
- Layers are always named and numbered within the data source
- Data sources can be directories containing files for the different layers or can be single files
- Data source can also be a relational database system or web service capable of providing simple features
vecdat = readOGR("PG:dbname=mydatabase", "thetable") # not run
This example would load the simple features and associated attribute columns from a PostGIS (http://www.postgis.org) table stored in a PostgreSQL database (http://www.postgresql.org).
Working with vector data
Vector data IO
plot(vecdat, col = "lightblue")
Working with vector data
Vector data IO
writeOGR does the opposite
dest = tempdir()
writeOGR(vecdat, dest, "newlayername", "ESRI Shapefile")
dir(dest)
unlink(dest) # remove the output
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
sp classes
class(vecdat) # sp class: help(package=sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
show(vecdat@data) # data slot contains attribute table
CONTINENT
0 Asia
1 North America
2 Europe
3 Africa
4 South America
5 Oceania
6 Australia
7 Antarctica
Working with vector data
sp classes
Classes
Spatially-referenced geometry | Geometries with attribute table
------------- | -------------
SpatialPoints | SpatialPointsDataFrame
SpatialLines | SpatialLinesDataFrame
SpatialPolygons | SpatialPolygonsDataFrame
The xDataFrame classes have a data frame attached. The rows of the data frame match the position in the list of geometries.
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Creating vector data
loc = rbind(c(30.25, 97.75), c(42.36, 71.06),
c(41.88, 87.63), c(37.78, 122.42))
loc = cbind(-loc[,2], loc[,1])
show(loc)
[,1] [,2]
[1,] -97.75 30.25
[2,] -71.06 42.36
[3,] -87.63 41.88
[4,] -122.42 37.78
Working with vector data
Creating vector data
spts = SpatialPoints(loc, CRS("+proj=longlat +ellps=WGS84"))
nam = subset(vecdat, CONTINENT=="North America")
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, xpd = NA, add = T) # plot continent boundaries
plot(spts, cex = 4, pch = 21, bg = "orange", add = T)
Working with vector data
Creating vector data
lut = data.frame(cities = c("Austin", "Boston",
"Chicago", "Los Angeles"))
spts.df = SpatialPointsDataFrame(loc, lut)
show(class(spts.df))
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
show(spts.df)
coordinates cities
1 (-97.75, 30.25) Austin
2 (-71.06, 42.36) Boston
3 (-87.63, 41.88) Chicago
4 (-122.42, 37.78) Los Angeles
Working with vector data
Creating vector data
loc = rbind(loc, loc[1,]) # 1st and last point same
lring = Polygon(loc, hole = FALSE)
poly = Polygons(list(lring), "1") # Single ring polygon
spoly = SpatialPolygons(list(poly), proj4 = CRS("+proj=longlat +ellps=WGS84"))
show(spoly)
An object of class "SpatialPolygons"
Slot "polygons":
[[1]]
An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] -96.76296 37.15718
Slot "area":
[1] 275.4836
Slot "hole":
[1] FALSE
Slot "ringDir":
[1] 1
Slot "coords":
[,1] [,2]
[1,] -97.75 30.25
[2,] -122.42 37.78
[3,] -87.63 41.88
[4,] -71.06 42.36
[5,] -97.75 30.25
Slot "plotOrder":
[1] 1
Slot "labpt":
[1] -96.76296 37.15718
Slot "ID":
[1] "1"
Slot "area":
[1] 275.4836
Slot "plotOrder":
[1] 1
Slot "bbox":
min max
x -122.42 -71.06
y 30.25 42.36
Slot "proj4string":
CRS arguments: +proj=longlat +ellps=WGS84
Working with vector data
Creating vector data
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(spoly, col = "lightblue", add = T)
plot(spts, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Plotting vector data
- Several different systems
- Base graphics
plot, points, lines, polygons- Grid graphics
- Low-level
- Base for lattice and ggplot2
- Lattice
xyplot, spplot- Ggplot2
qplot, ggplot- Lots of folks moving to this
Working with vector data
Plotting vector data
The plot command is specialized for sp classes and will generally do the right thing. A few switches:
type -- one of p, l, o, n (points, lines, points-lines, suppress output)add -- if TRUE then overlay on current plotcex -- symbol magnification factorlwd -- line widthcol -- the color to use (many ways to specify; see rgb)bg -- the background or fill colorxpd -- if set to NA will supress clipping at plot borderspch -- symbol for points (21 = filled circle)
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Reprojecting data
Working with vector data
Reprojecting data
p4s = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs" # US National Atlas
nam.laea = spTransform(nam, CRS(p4s))
spoly.laea = spTransform(spoly, CRS(p4s))
spts.laea = spTransform(spts, CRS(p4s))
show(spts.laea)
SpatialPoints:
coords.x1 coords.x2
[1,] 217894.6 -1632769.3
[2,] 2316750.9 124310.3
[3,] 1019668.8 -269814.5
[4,] -1945033.1 -538365.0
Coordinate Reference System (CRS) arguments: +proj=laea +lat_0=45
+lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
Working with vector data
Reprojecting data
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(spoly, col = "lightblue", add = T)
plot(spts, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
Reprojecting data
plot(spts.laea, type = "n") # set axis limits
plot(nam.laea, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(spoly.laea, col = "lightblue", add = T)
plot(spts.laea, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
Reprojecting data
- Based on PROJ4 software, see http://trac.osgeo.org/proj/
- Look up
proj4string at: http://spatialreference.org/
- Large libraries of "canned" definitions
- European Petroleum Survey Group (EPSG)
- ESRI
- International Astronomical Union (IAU)
- These define parameters, units, earth model, etc.
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
The rgeos library
- Wrapper for the Geometry Engine - Open Source
- C++ port of Java Topology Suite
- Used by PostGIS for geometry operations
- Also exposed by GDAL, but not wrapped in
rgdal (is in rgdal2)
Working with vector data
The rgeos library
Uses simple features
- POINT
- LINESTRING
- POLYGON
- MULTIPOINT
- MULTILINESTRING
- MULTIPOLYGON
- GEOMETRYCOLLECTION
stored in sp classes
Working with vector data
The rgeos library
Read in sp objecst using rgdal or other software or construct from WKT
pacman::p_load(rgeos)
p1 = readWKT("POINT(-97.75 30.25)", p4s = "+proj=longlat +ellps=WGS84")
show(p1)
SpatialPoints:
x y
1 -97.75 30.25
Coordinate Reference System (CRS) arguments: +proj=longlat
+ellps=WGS84
writeWKT(p1)
[1] "POINT (-97.7500000000000000 30.2500000000000000)"
Working with vector data
The rgeos library
p2 = readWKT("POINT(-71.06 42.36)", p4s = "+proj=longlat +ellps=WGS84")
gDistance(p1, p2) # Euclidian distance!
[1] 29.30884
pacman::p_load(geosphere)
distm(p1, p2) # in meters
| |
|-------:|
| 2731121|
gc = SpatialPoints(greatCircle(p1, p2), proj = CRS("+proj=longlat +ellps=WGS84"))
Working with vector data
The rgeos library
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(gc, lwd = 3, col = "lightblue", add = T)
plot(spts, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Geometry operators
Spatial sets
Working with vector data
Working with vector data
Geometry operators
Spatial sets
Working with vector data
Geometry operators
sp1 = readWKT("POLYGON ((10 10, 15 0, 25 0, 30 10, 25 20, 15 20, 10 10))")
sp2 = readWKT("POLYGON ((20 10, 30 0, 40 10, 30 20, 20 10))")
plot(gUnion(sp1, sp2), type = "n")
plot(sp1, col = "lightblue", add = T)
plot(sp2, col = "lightgreen", add = T)
plot(gIntersection(sp1, sp2), col = "lightyellow", add = T)
Working with vector data
Geometry operators
rel = gRelate(sp1, sp2); show(rel)
[1] "212101212"
mat = matrix(unlist(strsplit(rel, "")), 3)
typ = c("interior", "boundary", "exterior")
dimnames(mat) = list(typ, typ)
mode(mat) = "numeric"; show(mat)
interior boundary exterior
interior 2 1 2
boundary 1 0 1
exterior 2 1 2
Working with vector data
Geometry operators
Consider the following definition of Area/Area overlap:
OVERLAP
Interior
Boundary
Exterior
Interior
T
*
T
Boundary
*
*
*
Exterior
T
*
*
As a string: T\*T\*\*\*T\*\*
http://docs.geotools.org/stable/userguide/library/jts/dim9.html
Working with vector data
========================================================
Geometry operators
Relationship
Area/Area Pattern
“212101212”
Description
Disjoint
FF*FF****
false
x is not disjoint from y
Touches
FT*******
false
x does not just touch y
Touches
F***T****
false
x does not just touch y
Crosses
T*T***T**
true
x crosses y
Within
TF*F*****
false
x is not within y
Overlaps
T*T***T**
true
x overlaps y
Working with vector data
========================================================
Geometry operators
wzxhzdk:43
Working with vector data
========================================================
Geometry operators
wzxhzdk:44
Working with vector data
========================================================
Geometry operators
wzxhzdk:45
Working with vector data
========================================================
type: sub-section
- Simple features
- The `rgdal` package
- Vector data IO
- `sp` classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The `rgeos` library
- Geometry operators
- **Basic point processes**
- Interpolation and Kriging
Working with vector data
========================================================
Basic point processes
- Point process models focus on spatial pattern of locations
- Point patterns can be
- unmarked -- just the locations
- marked -- values associated with each point
- Main questions
- What is the intensity (points / area) of the pattern?
- Are the points related to each other in some way?
- Is the pattern the same everywhere or different in different places?
Working with vector data
========================================================
Basic point processes
- Libraries
- `spatial`
- `spatstat`
- `MarkedPointProcess`
- `splancs`
- many others
Working with vector data
========================================================
Basic point processes
wzxhzdk:46
wzxhzdk:47
wzxhzdk:48
wzxhzdk:49
Working with vector data
========================================================
Basic point processes
wzxhzdk:50
wzxhzdk:51
wzxhzdk:52
wzxhzdk:53
Working with vector data
========================================================
Basic point processes
wzxhzdk:54
Working with vector data
========================================================
Basic point processes
wzxhzdk:55
wzxhzdk:56
Working with vector data
========================================================
Basic point processes
wzxhzdk:57
wzxhzdk:58
wzxhzdk:59
wzxhzdk:60
wzxhzdk:61
wzxhzdk:62
Working with vector data
========================================================
Basic point processes
wzxhzdk:63
Working with vector data
========================================================
Basic point processes
wzxhzdk:64
Working with vector data
========================================================
Basic point processes
Working with vector data
========================================================
Basic point processes
Ripley's K $$\hat{K}(r) = \frac{a}{n(n-1)}\sum_i\sum_j I(d_{ij} < r) e_{ij}$$
Working with vector data
========================================================
Basic point processes
wzxhzdk:65
Working with vector data
========================================================
Basic point processes
Repeated Monte Carlo simulations of Poisson point process
wzxhzdk:66
wzxhzdk:67
Working with vector data
========================================================
Basic point processes
wzxhzdk:68
Working with vector data
========================================================
Basic point processes
wzxhzdk:69
wzxhzdk:70
Working with vector data
========================================================
Basic point processes
wzxhzdk:71
Working with vector data
========================================================
Basic point processes
wzxhzdk:72
Working with vector data
========================================================
type: sub-section
- Simple features
- The `rgdal` package
- Vector data IO
- `sp` classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The `rgeos` library
- Geometry operators
- Basic point processes
- **Interpolation and Kriging**
Working with vector data
========================================================
Interpolation and Kriging
Thin-plate spline:
$$ y_i = f(x_i) + \epsilon_i $$
where $x_i$ are coordinates, $y_i$ are measurments and $f$ is a smooth function found by minimizing:
$$||\mathbf{y}-\mathbf{f}||^2 + \lambda \int \mathbf{f}''(x)^2 dx$$
where $\lambda$ is a smoothing parameter. The resulting $f$ smoothly interpolates the data. It is an esimate of the locally conditioned mean of $\mathbf{y}$.
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:73
wzxhzdk:74
wzxhzdk:75
wzxhzdk:76
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:77
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:78
wzxhzdk:79
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:80
Working with vector data
========================================================
Interpolation and Kriging
Regression Kriging model
$$ Y_k = P(x_k) + Z(x_k) + \epsilon_k $$
where $P(x)$ is a low-order polynomial trend surface and $Z(x)$ is a stationary Gaussian processes with spatial covariance $\Sigma(||x_i-x_j||)$. A thin-plate spline is a special case of Kriging -- the `fields` package function `Tps` actually calls `Krig` internally.
Many more packages. See `gstat`, `geoR` and `geoRglm` among others.
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:81
wzxhzdk:82
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:83
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:84
Working with vector data
========================================================
type: sub-section
- Simple features
- The `rgdal` package
- Vector data IO
- `sp` classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The `rgeos` library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
thk686/keitt.ssi.2019 documentation built on June 28, 2019, 1:28 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Geospatial Data Analysis in R
author: Timothy H. Keitt date: May 12, 2014 width: 1440 height: 900
Working with vector data
type: section
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Spatial layers
GIS building blocks
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
OGC Simple Feature Hierarchy
Working with vector data
OGC Simple Feature Hierarchy
Working with vector data
Complex versus simple features
Most spatial operators in R require simple features
Working with vector data
Complex versus simple features
Most spatial operators in R require simple features
Working with vector data
OGC Simple Features Well-Known-Text
Working with vector data
OGC Simple Features Well-Known-Text
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
RGDAL package
- OGR component of GDAL handles vector data (http://www.gdal.org/)
- Data source -> data layer -> simple features
- Written in C++
- Like GDAL, large collection of runtime drivers
- Encapsulated in
rgdaland experimentalrgdal2packages - With
rgdalpackage, main use is to read data intospclasses
Working with vector data
RGDAL package
pacman::p_load(rgdal)
as.character(ogrDrivers()[[1]])
[1] "AeronavFAA" "AmigoCloud" "ARCGEN" "AVCBin"
[5] "AVCE00" "BNA" "Carto" "Cloudant"
[9] "CouchDB" "CSV" "CSW" "DGN"
[13] "DXF" "EDIGEO" "ElasticSearch" "ESRI Shapefile"
[17] "Geoconcept" "GeoJSON" "GeoRSS" "GFT"
[21] "GML" "GPKG" "GPSBabel" "GPSTrackMaker"
[25] "GPX" "HTF" "HTTP" "Idrisi"
[29] "JML" "KML" "MapInfo File" "Memory"
[33] "netCDF" "ODS" "OGR_GMT" "OGR_PDS"
[37] "OGR_SDTS" "OGR_VRT" "OpenAir" "OpenFileGDB"
[41] "OSM" "PCIDSK" "PDF" "PGDUMP"
[45] "PLSCENES" "REC" "S57" "SEGUKOOA"
[49] "SEGY" "Selafin" "SQLite" "SUA"
[53] "SVG" "SXF" "TIGER" "UK .NTF"
[57] "VDV" "VFK" "WAsP" "WFS"
[61] "XLSX" "XPlane"
Uniform access to many data sources
Uniform access to many data sources
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Vector data IO
- main functions
- readOGR
- writeOGR
vecdat = readOGR("example-data/continents", "continent")
OGR data source with driver: ESRI Shapefile
Source: "/Users/tkeitt/Dropbox/R/keitt.ssi.2019/inst/materials/lectures/example-data/continents", layer: "continent"
with 8 features
It has 1 fields
First argument (dsn) is the data source. The second argument is the layer within the data source.
Working with vector data
Vector data IO
An ESRI shapefile is often several files stored together in a directory. The directory is the data source and the .shp file is the layer. Other files hold spatial reference system and tabular field data.
dir("example-data") # see the data source
[1] "continents" "NCEAS sample" "ozone.gml" "ozone.xsd"
[5] "rainforest" "rainforest.tif"
dir("example-data/continents") # see the data layers
[1] "continent.dbf" "continent.prj" "continent.shp" "continent.shx"
Other formats are single file, in which case the data source is the file name and the layer is named within the file, typically the same as the file name without the extension.
Working with vector data
Vector data IO
- Layers are always named and numbered within the data source
- Data sources can be directories containing files for the different layers or can be single files
- Data source can also be a relational database system or web service capable of providing simple features
vecdat = readOGR("PG:dbname=mydatabase", "thetable") # not run
This example would load the simple features and associated attribute columns from a PostGIS (http://www.postgis.org) table stored in a PostgreSQL database (http://www.postgresql.org).
Working with vector data
Vector data IO
plot(vecdat, col = "lightblue")
Working with vector data
Vector data IO
writeOGR does the opposite
dest = tempdir()
writeOGR(vecdat, dest, "newlayername", "ESRI Shapefile")
dir(dest)
unlink(dest) # remove the output
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
sp classes
class(vecdat) # sp class: help(package=sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
show(vecdat@data) # data slot contains attribute table
CONTINENT
0 Asia
1 North America
2 Europe
3 Africa
4 South America
5 Oceania
6 Australia
7 Antarctica
Working with vector data
sp classes
Classes
Spatially-referenced geometry | Geometries with attribute table ------------- | ------------- SpatialPoints | SpatialPointsDataFrame SpatialLines | SpatialLinesDataFrame SpatialPolygons | SpatialPolygonsDataFrame
The xDataFrame classes have a data frame attached. The rows of the data frame match the position in the list of geometries.
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Creating vector data
loc = rbind(c(30.25, 97.75), c(42.36, 71.06),
c(41.88, 87.63), c(37.78, 122.42))
loc = cbind(-loc[,2], loc[,1])
show(loc)
[,1] [,2]
[1,] -97.75 30.25
[2,] -71.06 42.36
[3,] -87.63 41.88
[4,] -122.42 37.78
Working with vector data
Creating vector data
spts = SpatialPoints(loc, CRS("+proj=longlat +ellps=WGS84"))
nam = subset(vecdat, CONTINENT=="North America")
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, xpd = NA, add = T) # plot continent boundaries
plot(spts, cex = 4, pch = 21, bg = "orange", add = T)
Working with vector data
Creating vector data
lut = data.frame(cities = c("Austin", "Boston",
"Chicago", "Los Angeles"))
spts.df = SpatialPointsDataFrame(loc, lut)
show(class(spts.df))
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
show(spts.df)
coordinates cities
1 (-97.75, 30.25) Austin
2 (-71.06, 42.36) Boston
3 (-87.63, 41.88) Chicago
4 (-122.42, 37.78) Los Angeles
Working with vector data
Creating vector data
loc = rbind(loc, loc[1,]) # 1st and last point same
lring = Polygon(loc, hole = FALSE)
poly = Polygons(list(lring), "1") # Single ring polygon
spoly = SpatialPolygons(list(poly), proj4 = CRS("+proj=longlat +ellps=WGS84"))
show(spoly)
An object of class "SpatialPolygons"
Slot "polygons":
[[1]]
An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] -96.76296 37.15718
Slot "area":
[1] 275.4836
Slot "hole":
[1] FALSE
Slot "ringDir":
[1] 1
Slot "coords":
[,1] [,2]
[1,] -97.75 30.25
[2,] -122.42 37.78
[3,] -87.63 41.88
[4,] -71.06 42.36
[5,] -97.75 30.25
Slot "plotOrder":
[1] 1
Slot "labpt":
[1] -96.76296 37.15718
Slot "ID":
[1] "1"
Slot "area":
[1] 275.4836
Slot "plotOrder":
[1] 1
Slot "bbox":
min max
x -122.42 -71.06
y 30.25 42.36
Slot "proj4string":
CRS arguments: +proj=longlat +ellps=WGS84
Working with vector data
Creating vector data
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(spoly, col = "lightblue", add = T)
plot(spts, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Plotting vector data
- Several different systems
- Base graphics
plot,points,lines,polygons- Grid graphics
- Low-level
- Base for lattice and ggplot2
- Lattice
xyplot,spplot- Ggplot2
qplot,ggplot- Lots of folks moving to this
Working with vector data
Plotting vector data
The plot command is specialized for sp classes and will generally do the right thing. A few switches:
type-- one ofp, l, o, n(points, lines, points-lines, suppress output)add-- ifTRUEthen overlay on current plotcex-- symbol magnification factorlwd-- line widthcol-- the color to use (many ways to specify; seergb)bg-- the background or fill colorxpd-- if set toNAwill supress clipping at plot borderspch-- symbol for points (21 = filled circle)
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Reprojecting data
Working with vector data
Reprojecting data
p4s = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs" # US National Atlas
nam.laea = spTransform(nam, CRS(p4s))
spoly.laea = spTransform(spoly, CRS(p4s))
spts.laea = spTransform(spts, CRS(p4s))
show(spts.laea)
SpatialPoints:
coords.x1 coords.x2
[1,] 217894.6 -1632769.3
[2,] 2316750.9 124310.3
[3,] 1019668.8 -269814.5
[4,] -1945033.1 -538365.0
Coordinate Reference System (CRS) arguments: +proj=laea +lat_0=45
+lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
Working with vector data
Reprojecting data
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(spoly, col = "lightblue", add = T)
plot(spts, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
Reprojecting data
plot(spts.laea, type = "n") # set axis limits
plot(nam.laea, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(spoly.laea, col = "lightblue", add = T)
plot(spts.laea, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
Reprojecting data
- Based on PROJ4 software, see http://trac.osgeo.org/proj/
- Look up
proj4stringat: http://spatialreference.org/ - Large libraries of "canned" definitions
- European Petroleum Survey Group (EPSG)
- ESRI
- International Astronomical Union (IAU)
- These define parameters, units, earth model, etc.
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
The rgeos library
- Wrapper for the Geometry Engine - Open Source
- C++ port of Java Topology Suite
- Used by PostGIS for geometry operations
- Also exposed by GDAL, but not wrapped in
rgdal(is inrgdal2)
Working with vector data
The rgeos library
Uses simple features
- POINT
- LINESTRING
- POLYGON
- MULTIPOINT
- MULTILINESTRING
- MULTIPOLYGON
- GEOMETRYCOLLECTION
stored in sp classes
Working with vector data
The rgeos library
Read in sp objecst using rgdal or other software or construct from WKT
pacman::p_load(rgeos)
p1 = readWKT("POINT(-97.75 30.25)", p4s = "+proj=longlat +ellps=WGS84")
show(p1)
SpatialPoints:
x y
1 -97.75 30.25
Coordinate Reference System (CRS) arguments: +proj=longlat
+ellps=WGS84
writeWKT(p1)
[1] "POINT (-97.7500000000000000 30.2500000000000000)"
Working with vector data
The rgeos library
p2 = readWKT("POINT(-71.06 42.36)", p4s = "+proj=longlat +ellps=WGS84")
gDistance(p1, p2) # Euclidian distance!
[1] 29.30884
pacman::p_load(geosphere)
distm(p1, p2) # in meters
| | |-------:| | 2731121|
gc = SpatialPoints(greatCircle(p1, p2), proj = CRS("+proj=longlat +ellps=WGS84"))
Working with vector data
The rgeos library
plot(spts, type = "n") # set axis limits
plot(nam, lwd = 2, col = "navajowhite", xpd = NA, add = T)
plot(gc, lwd = 3, col = "lightblue", add = T)
plot(spts, cex = 4, pch = 21, bg = "darkgreen", add = T)
Working with vector data
type: sub-section
- Simple features
- The rgdal package
- Vector data IO
- sp classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The rgeos library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
Working with vector data
Geometry operators
Spatial sets
Working with vector data
Working with vector data
Geometry operators
Spatial sets
Working with vector data
Geometry operators
sp1 = readWKT("POLYGON ((10 10, 15 0, 25 0, 30 10, 25 20, 15 20, 10 10))")
sp2 = readWKT("POLYGON ((20 10, 30 0, 40 10, 30 20, 20 10))")
plot(gUnion(sp1, sp2), type = "n")
plot(sp1, col = "lightblue", add = T)
plot(sp2, col = "lightgreen", add = T)
plot(gIntersection(sp1, sp2), col = "lightyellow", add = T)
Working with vector data
Geometry operators
rel = gRelate(sp1, sp2); show(rel)
[1] "212101212"
mat = matrix(unlist(strsplit(rel, "")), 3)
typ = c("interior", "boundary", "exterior")
dimnames(mat) = list(typ, typ)
mode(mat) = "numeric"; show(mat)
interior boundary exterior
interior 2 1 2
boundary 1 0 1
exterior 2 1 2
Working with vector data
Geometry operators
Consider the following definition of Area/Area overlap:
OVERLAP Interior Boundary Exterior Interior T * T Boundary * * * Exterior T * *As a string: T\*T\*\*\*T\*\* http://docs.geotools.org/stable/userguide/library/jts/dim9.html Working with vector data ======================================================== Geometry operators
Relationship Area/Area Pattern “212101212” Description Disjoint FF*FF**** false x is not disjoint from y Touches FT******* false x does not just touch y Touches F***T**** false x does not just touch y Crosses T*T***T** true x crosses y Within TF*F***** false x is not within y Overlaps T*T***T** true x overlaps y  Working with vector data ======================================================== Geometry operators wzxhzdk:43  Working with vector data ======================================================== Geometry operators wzxhzdk:44  Working with vector data ======================================================== Geometry operators wzxhzdk:45  Working with vector data ======================================================== type: sub-section - Simple features - The `rgdal` package - Vector data IO - `sp` classes - Creating vector data - Plotting vector data - Reprojecting data - The `rgeos` library - Geometry operators - **Basic point processes** - Interpolation and Kriging Working with vector data ======================================================== Basic point processes - Point process models focus on spatial pattern of locations - Point patterns can be - unmarked -- just the locations - marked -- values associated with each point - Main questions - What is the intensity (points / area) of the pattern? - Are the points related to each other in some way? - Is the pattern the same everywhere or different in different places? Working with vector data ======================================================== Basic point processes - Libraries - `spatial` - `spatstat` - `MarkedPointProcess` - `splancs` - many others Working with vector data ======================================================== Basic point processes wzxhzdk:46 wzxhzdk:47 wzxhzdk:48 wzxhzdk:49 Working with vector data ======================================================== Basic point processes wzxhzdk:50 wzxhzdk:51 wzxhzdk:52 wzxhzdk:53 Working with vector data ======================================================== Basic point processes wzxhzdk:54  Working with vector data ======================================================== Basic point processes wzxhzdk:55 wzxhzdk:56 Working with vector data ======================================================== Basic point processes wzxhzdk:57 wzxhzdk:58 wzxhzdk:59 wzxhzdk:60 wzxhzdk:61 wzxhzdk:62 Working with vector data ======================================================== Basic point processes wzxhzdk:63  Working with vector data ======================================================== Basic point processes wzxhzdk:64 Working with vector data ======================================================== Basic point processes  Working with vector data ======================================================== Basic point processes Ripley's K $$\hat{K}(r) = \frac{a}{n(n-1)}\sum_i\sum_j I(d_{ij} < r) e_{ij}$$
Working with vector data
========================================================
Basic point processes
wzxhzdk:65
Working with vector data
========================================================
Basic point processes
Repeated Monte Carlo simulations of Poisson point process
wzxhzdk:66
wzxhzdk:67
Working with vector data
========================================================
Basic point processes
wzxhzdk:68
Working with vector data
========================================================
Basic point processes
wzxhzdk:69
wzxhzdk:70
Working with vector data
========================================================
Basic point processes
wzxhzdk:71
Working with vector data
========================================================
Basic point processes
wzxhzdk:72
Working with vector data
========================================================
type: sub-section
- Simple features
- The `rgdal` package
- Vector data IO
- `sp` classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The `rgeos` library
- Geometry operators
- Basic point processes
- **Interpolation and Kriging**
Working with vector data
========================================================
Interpolation and Kriging
Thin-plate spline:
$$ y_i = f(x_i) + \epsilon_i $$
where $x_i$ are coordinates, $y_i$ are measurments and $f$ is a smooth function found by minimizing:
$$||\mathbf{y}-\mathbf{f}||^2 + \lambda \int \mathbf{f}''(x)^2 dx$$
where $\lambda$ is a smoothing parameter. The resulting $f$ smoothly interpolates the data. It is an esimate of the locally conditioned mean of $\mathbf{y}$.
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:73
wzxhzdk:74
wzxhzdk:75
wzxhzdk:76
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:77
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:78
wzxhzdk:79
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:80
Working with vector data
========================================================
Interpolation and Kriging
Regression Kriging model
$$ Y_k = P(x_k) + Z(x_k) + \epsilon_k $$
where $P(x)$ is a low-order polynomial trend surface and $Z(x)$ is a stationary Gaussian processes with spatial covariance $\Sigma(||x_i-x_j||)$. A thin-plate spline is a special case of Kriging -- the `fields` package function `Tps` actually calls `Krig` internally.
Many more packages. See `gstat`, `geoR` and `geoRglm` among others.
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:81
wzxhzdk:82
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:83
Working with vector data
========================================================
Interpolation and Kriging
wzxhzdk:84
Working with vector data
========================================================
type: sub-section
- Simple features
- The `rgdal` package
- Vector data IO
- `sp` classes
- Creating vector data
- Plotting vector data
- Reprojecting data
- The `rgeos` library
- Geometry operators
- Basic point processes
- Interpolation and Kriging
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.