What is a GIS?
Software and file types
Projections and coordinate reference systems
Spatially referenced data in R
Sourcing
Manipulating
Visualizing
" system designed to capture, store, manipulate, analyze, manage, and present all types of spatial or geographical data"
-- wiki
Term coined in 1968
-Roger Tomlinson, Canadian regional planner
Relational database where index field is spatiotemporal reference
Generally, baselines are from topos, enhanced with satellite imagery / aerial data
Vector
Raster
3D color images
Source: bgis.sanbi.org
ESRI's file format
Source: bgis.sanbi.org
ArcGIS (ESRI)
Google Maps / Earth (KML)
GDAL/OGR (Geographic data abstraction library)
GEOS (topological projections)
KML = "keynote markup language"
Spatial vsn of XML
Systematic transformation of locations on a sphere to locations on a plane
Developable surfaces
For sphere, some metric always gets distorted
Source: Geographic Information System Basics, Schmitz (2012)
1. Project to developable surface
2. Unroll surface into plane
Flattened to plane
Developable surfaces
Source: http://www.progonos.com
What metric to distort?
Check out all these projections!
Examples
Elements
rgdal ------------------------ brings GDAL/OGR functionality to R
spatial ---------------------- defines key spatial object classes
sp ---------------------------- transformation and drilling functionality
maptools ------------------ gets a lot of use, especially readShapePoly
rgeos ----------------------- interface engine connecting R with GEOS
(C API with lots of nice topological functionality)
(Spatial Points (Polys, Etc) Data Frames)
1. Read in csv with coordinates stored in "Long" and "Lat"
MyData <- read.csv("<Path to data.csv>", header = T)
class(MyData) == data.frame
2. Use coordinates() function to specify coordinates
coordinates(MyData) = ~ Long + Lat
class(MyData) == SPDF
3. Specify coordinate reference system
latlong <- "+init=epsg:4326"
proj4string(MyData) <- CRS(latlong)
Political boundaries
Topography
Lots of map tiles available through leaflet
Enter google maps
Navigate to location of interest
Get coordinates:
require(RgoogleMaps)
center.in <- c(45.358, -111.403)
map <- GetMap(center = center.in, zoom = 13, maptype = "terrain")
maptools
rgdal
with maptools
GoogleEarthData <- getKMLcoordinates(kmlfile = "./Data/MyPoints.kml")
Points come in a list, so
GoogleEarthData <- do.call("rbind", GoogleEarthData)
Check current projection with
proj4string(MySPDF)
spTransform() in sp to switch to a new projection
MySPDFReprojected <- spTransform(MySPDF, CRS(latlong))
(recall: we defined latlong earlier to be an object specifying a particular CRS)
Load RgoogleMaps
require(RgoogleMaps)
Project data to decimal-degrees if it's not already there
MySPDFLatLong <- spTransform(MySPDF, CRS(latlong))
Recenter lat/longs relative to google map tile
MySPDFGoogle <- LatLon2XY.centered(MyMapTile,
MySPDFLatLong$Latitude,
MySPDFLatLong$Longitude,
zoom = GoogleZoom)
Goal:
Align two layers & extracts values from one layer that occur at spatial locations specified in the other
over() in sp
MyPoints <- readOGR(dsn = "./PathToPoints", layer = "MyPointsLayer")
MyPolys <- readOGR(dsn = "./PathToPolys", layer = "MyPolysLayer")
PolyValsAtPoints <- (MyPoints %over% MyPolys)$Value
Scenario:
Shapefile of polygons; extract overlap area for each pair
gIntersection() in rgeos
MyPolys <- readOGR(dsn = "./PathToPolys", layer = "MyPolysLayer")
# Extract geometries of the two focal polygons
Poly1 <- MyPolys$Polygons[[1]] # gets 1st polygon from MyPolys
Poly2 <- MyPolys$Polygons[[2]] # gets 2nd polygon from MyPolys
require(rgeos)
Poly1Poly2Overlap <- gIntersection(Poly1, Poly2, byid = F)
1. Import map tile from google maps
require(RgoogleMaps)
center.in <- c(45.358, -111.403)
MyMap <- GetMap(center = center.in,
zoom = googlezoom,
maptype = "terrain")
2. Project points and store in dataframe
MyPoints <- readOGR("./PathToPoints", layer = "MyPoints")
MyPointsGoogle <- LatLon2XY.centered(MyMap,
MyPoints$Lat, MyPoints$Long, zoom = googlezoom)
MyPointsGoogleData <- data.frame(MyPointsGoogle)
3. Assemble plot
PlotOnStaticMap(MyMap)
points(MyPointsGoogleData$newX, MyPointsGoogleData$newY,
pch = 16, col = "blue")