Shapefiles and Rasters
Naomi Schwartz and Dan McGlinn
# install.packages(c("maps","sf", "raster"))
library(maps) # convinient pkg for maps of the world
library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE# First, make some maps. The "maps" package has databases="world","usa","state",
# or "county". (There are also a few for foreign countries like France and
# Italy.) In each database, you can plot any region or group of regions.
world <- map(database="world")
names(world)## [1] "x" "y" "range" "names"head(world$names)## [1] "Aruba" "Afghanistan" "Angola" "Angola:Cabinda"
## [5] "Anguilla" "Albania"# names in world database can all be plot as a region. For example:
map(database="world", regions=c("Cambodia","Thailand","Vietnam","Laos"))
map(database="state",regions=c("virginia","north carolina","south carolina"))
# you can add cities using map.cities command.
nc_county <- map(database="county",regions="north carolina")
head(nc_county$names)## [1] "north carolina,alamance" "north carolina,alexander"
## [3] "north carolina,alleghany" "north carolina,anson"
## [5] "north carolina,ashe" "north carolina,avery"triangle <- map(database="county",
regions=c("north carolina,orange","north carolina,wake","north carolina,durham"),
fill=T, plot=F)
triangle$names## [1] "north carolina,durham" "north carolina,orange" "north carolina,wake"# We can turn a map into a spatial polygons object and fill it with data,
# turning it into a spatial polygon dataframe.
tri_sf <- st_as_sf(triangle)
plot(tri_sf, axes=T)
# add data to make spatial polygons data frame
population <- c(266132, 132272, 892409)
tri_sf$pop <- population
plot(tri_sf['pop'])
# If your data are in different projections, you need to change the projection
# so that they are all in
# the same coordinate reference system.
world <- map(database="world", fill=T, plot=F)
world_longlat <- st_as_sf(world)
# use fucntion spTransform to transform to mercator
world_merc <- st_transform(world_longlat, crs=st_crs("+proj=merc"))
# Lambert Azimuthal Equal Area
world_laea <- st_transform(world_longlat, crs=st_crs("+proj=laea"))
# sinusoidal
world_sinusoidal <- st_transform(world_longlat,
crs=st_crs("+proj=sinu"))
# plot the four together to see the difference. Mercator projection distorts
# area far from the equator. LAEA is accurately represents area but not angles.
# Sinusoidal is equal area and conserves distances along parallels.
par(mfrow=c(1,1))
plot(world_longlat)
plot(world_merc)
plot(world_laea)
plot(world_sinusoidal)
# A four page cheat sheet about coordinate reference systems (CRSs), including
# projections, datums, and coordinate systems, and the use of these in R.
#https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf
# A fun projection link
# http://xkcd.com/977/
# OK, some slightly more complicated stuff. Data from MODIS fire detections in
# the US in 2010.
# metadata for MODIS fire data available at:
# http://firemapper.sc.egov.usda.gov/data/fireptdata/modisfire_2010_conus.htm
fire2010 <- sf::st_read(dsn = './data/MODISfire2010',
layer = 'modis_fire_2010_365_conus')## Reading layer `modis_fire_2010_365_conus' from data source
## `/home/mcglinndj/quant_methods/data/MODISfire2010' using driver `ESRI Shapefile'
## Simple feature collection with 131586 features and 16 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -124.883 ymin: 25.165 xmax: -67.379 ymax: 49.454
## Geodetic CRS: NAD83class(fire2010)## [1] "sf" "data.frame"names(fire2010) ## [1] "AREA" "PERIMETER" "FIRE_" "FIRE_ID" "LAT" "LONG"
## [7] "DATE" "JULIAN" "GMT" "TEMP" "SPIX" "TPIX"
## [13] "SRC" "SAT_SRC" "CONF" "FRP" "geometry"st_crs(fire2010) # retrieves coordinate reference system## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## 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",4269]]dim(fire2010)## [1] 131586 17# make a map of Julian day of each fire
plot(fire2010["JULIAN"], cex = 0.25)
# get summary statistics for each attribute
summary(fire2010)## AREA PERIMETER FIRE_ FIRE_ID LAT
## Min. :0 Min. :0 Min. : 1 Min. : 63 Min. :25.16
## 1st Qu.:0 1st Qu.:0 1st Qu.: 32897 1st Qu.: 76054 1st Qu.:32.02
## Median :0 Median :0 Median : 65794 Median :302178 Median :34.91
## Mean :0 Mean :0 Mean : 65794 Mean :249912 Mean :36.01
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.: 98690 3rd Qu.:391697 3rd Qu.:38.70
## Max. :0 Max. :0 Max. :131586 Max. :482122 Max. :49.45
## LONG DATE JULIAN GMT
## Min. :-124.88 Min. :2010-01-01 Min. : 1.0 Min. : 233
## 1st Qu.:-100.07 1st Qu.:2010-04-16 1st Qu.:106.0 1st Qu.:1805
## Median : -92.73 Median :2010-08-31 Median :243.0 Median :1858
## Mean : -95.71 Mean :2010-07-27 Mean :208.4 Mean :1803
## 3rd Qu.: -87.06 3rd Qu.:2010-10-08 3rd Qu.:281.0 3rd Qu.:1931
## Max. : -67.38 Max. :2010-12-31 Max. :365.0 Max. :2210
## TEMP SPIX TPIX SRC
## Min. :305.0 Min. :1.000 Min. :1.000 Length:131586
## 1st Qu.:315.0 1st Qu.:1.100 1st Qu.:1.000 Class :character
## Median :320.9 Median :1.300 Median :1.100 Mode :character
## Mean :326.1 Mean :1.564 Mean :1.193
## 3rd Qu.:330.9 3rd Qu.:1.700 3rd Qu.:1.300
## Max. :502.6 Max. :4.800 Max. :2.000
## SAT_SRC CONF FRP geometry
## Length:131586 Min. : 0.00 Min. :-9999.00 POINT :131586
## Class :character 1st Qu.: 56.00 1st Qu.: 16.60 epsg:4269 : 0
## Mode :character Median : 70.00 Median : 29.60 +proj=long...: 0
## Mean : 68.52 Mean : 44.01
## 3rd Qu.: 83.00 3rd Qu.: 60.10
## Max. :100.00 Max. : 5317.30# subsetting: select fires that occurred in the last 6 months (days 182 to 365)
# of 2010
fire <- subset(fire2010, JULIAN > 182)
class(fire)## [1] "sf" "data.frame"names(fire)## [1] "AREA" "PERIMETER" "FIRE_" "FIRE_ID" "LAT" "LONG"
## [7] "DATE" "JULIAN" "GMT" "TEMP" "SPIX" "TPIX"
## [13] "SRC" "SAT_SRC" "CONF" "FRP" "geometry"dim(fire)## [1] 82802 17summary(fire)## AREA PERIMETER FIRE_ FIRE_ID LAT
## Min. :0 Min. :0 Min. : 4 Min. :190935 Min. :25.36
## 1st Qu.:0 1st Qu.:0 1st Qu.: 32126 1st Qu.:309123 1st Qu.:32.31
## Median :0 Median :0 Median : 65616 Median :356566 Median :34.94
## Mean :0 Mean :0 Mean : 64651 Mean :359921 Mean :36.19
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.: 95393 3rd Qu.:418602 3rd Qu.:38.93
## Max. :0 Max. :0 Max. :131580 Max. :482122 Max. :49.45
## LONG DATE JULIAN GMT
## Min. :-124.88 Min. :2010-07-02 Min. :183 Min. : 233
## 1st Qu.:-110.42 1st Qu.:2010-09-06 1st Qu.:249 1st Qu.:1700
## Median : -92.50 Median :2010-09-30 Median :273 Median :1840
## Mean : -97.42 Mean :2010-10-01 Mean :274 Mean :1752
## 3rd Qu.: -89.33 3rd Qu.:2010-10-23 3rd Qu.:296 3rd Qu.:1923
## Max. : -68.04 Max. :2010-12-31 Max. :365 Max. :2210
## TEMP SPIX TPIX SRC
## Min. :305.0 Min. :1.000 Min. :1.000 Length:82802
## 1st Qu.:315.2 1st Qu.:1.100 1st Qu.:1.000 Class :character
## Median :321.4 Median :1.300 Median :1.100 Mode :character
## Mean :326.7 Mean :1.571 Mean :1.195
## 3rd Qu.:331.6 3rd Qu.:1.700 3rd Qu.:1.300
## Max. :502.6 Max. :4.800 Max. :2.000
## SAT_SRC CONF FRP geometry
## Length:82802 Min. : 0.00 Min. : 0.00 POINT :82802
## Class :character 1st Qu.: 57.00 1st Qu.: 17.10 epsg:4269 : 0
## Mode :character Median : 71.00 Median : 30.90 +proj=long...: 0
## Mean : 69.47 Mean : 62.36
## 3rd Qu.: 84.00 3rd Qu.: 62.60
## Max. :100.00 Max. :5317.30plot(fire['JULIAN'], cex = 0.25)
# select the fires in NC. First we'll identify which state each fire occurred
# in, then subset the ones from NC.
USA <- map(database='state', plot=F, fill=T)
names(USA)## [1] "x" "y" "range" "names"USA$names## [1] "alabama" "arizona"
## [3] "arkansas" "california"
## [5] "colorado" "connecticut"
## [7] "delaware" "district of columbia"
## [9] "florida" "georgia"
## [11] "idaho" "illinois"
## [13] "indiana" "iowa"
## [15] "kansas" "kentucky"
## [17] "louisiana" "maine"
## [19] "maryland" "massachusetts:martha's vineyard"
## [21] "massachusetts:main" "massachusetts:nantucket"
## [23] "michigan:north" "michigan:south"
## [25] "minnesota" "mississippi"
## [27] "missouri" "montana"
## [29] "nebraska" "nevada"
## [31] "new hampshire" "new jersey"
## [33] "new mexico" "new york:manhattan"
## [35] "new york:main" "new york:staten island"
## [37] "new york:long island" "north carolina:knotts"
## [39] "north carolina:main" "north carolina:spit"
## [41] "north dakota" "ohio"
## [43] "oklahoma" "oregon"
## [45] "pennsylvania" "rhode island"
## [47] "south carolina" "south dakota"
## [49] "tennessee" "texas"
## [51] "utah" "vermont"
## [53] "virginia:chesapeake" "virginia:chincoteague"
## [55] "virginia:main" "washington:san juan island"
## [57] "washington:lopez island" "washington:orcas island"
## [59] "washington:whidbey island" "washington:main"
## [61] "west virginia" "wisconsin"
## [63] "wyoming"USA_sf <- st_as_sf(USA)
# make sure projection of USA polygon is same as fire points
USA_sf <- st_transform(USA_sf, st_crs(fire))
# overlay performs a "point in polygon" operation--meaning that it will return
# us a vector giving the index of which polygon in USA_sp each point in firehot
# is. We then index that to names(USA_sp) to get the name of the state for that
# index.
sf_use_s2(FALSE)## Spherical geometry (s2) switched offfire_state <- st_intersection(fire, USA_sf)## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# Then, subset the ones from NC. We have to use grep because there are actually
# three polygons for NC.
firenc <- fire_state[grep("north carolina", fire_state$ID), ]
dim(firenc)## [1] 1188 18# make a map of fires by temperature.
addLegendToSFPlot <- function(values = c(0, 1), labels = c("Low", "High"),
palette = c("blue", "red"), ...){
# Get the axis limits and calculate size
axisLimits <- par()$usr
xLength <- axisLimits[2] - axisLimits[1]
yLength <- axisLimits[4] - axisLimits[3]
# Define the colour palette
colourPalette <- leaflet::colorNumeric(palette, range(values))
# Add the legend
plotrix::color.legend(xl=axisLimits[2] - 0.1*xLength, xr=axisLimits[2],
yb=axisLimits[3], yt=axisLimits[3] + 0.1 * yLength,
legend = labels, rect.col = colourPalette(values),
gradient="y", ...)
}
#First make an NC spatial polygon to put
# around the fire points.
NC <- map(database='state', regions='north carolina', fill=T,
plot=F)
NC_st <- st_as_sf(NC)
NC_st <- st_transform(NC_st, crs=st_crs(firenc))
cuts <- cut(firenc$JULIAN, 10)
colors <- heat.colors(10)[as.numeric(cuts)]
plot(st_geometry(NC_st))
plot(st_geometry(firenc['JULIAN']), add = TRUE,
col = colors, pch = 19, cex = 0.5)
addLegendToSFPlot(values = seq(from = 183, to = 363, length.out = 10),
labels = c("Low", "", "", "", "Medium", "", "", "", "", "High"),
palette = heat.colors(10))
# cannot figure out how to combine in ggplot nicely.
#ggplot(NC_st) + geom_sf()
#ggplot(firenc['JULIAN']) + geom_sf()
#geom_point(st_geometry(firenc['JULIAN']))
# for the last step, let's export this as a KML (readable by google earth) using
# write OGR and plot the locations of fires in NC in google earth. The first
# argument is the object we want to export, the second is the filename (by
# default it will go in our working directory), the layer we want to export, and
# the file format.
write_sf(firenc, "firenctemp.kml", "TEMP", driver="kml")
## Rasters ----------------------------------------
## download Rdata file here:
## https://www.dropbox.com/s/gafxazc9575nf3j/bioclim_10m.Rdata?dl=0
library(raster)## Loading required package: spload('./data/bioclim_10m.Rdata')
plot(bioStack, "mat")
?extract## Help on topic 'extract' was found in the following packages:
##
## Package Library
## magrittr /home/mcglinndj/R/x86_64-pc-linux-gnu-library/4.3
## raster /home/mcglinndj/R/x86_64-pc-linux-gnu-library/4.3
## terra /home/mcglinndj/R/x86_64-pc-linux-gnu-library/4.3
##
##
## Using the first match ...fire_climate <- extract(bioStack, firenc)
class(fire_climate)## [1] "matrix" "array"head(fire_climate)## mat mdr iso tseas tmax tmin tar twetq tdryq twarmq tcoldq ap pwet
## [1,] 14.1 13.5 3.8 787.9 31.7 -3.3 35.0 24.1 5.4 24.1 3.7 1112 111
## [2,] 11.3 13.2 3.9 742.2 27.8 -5.6 33.4 18.9 1.4 20.6 1.4 1245 122
## [3,] 11.3 13.2 3.9 742.2 27.8 -5.6 33.4 18.9 1.4 20.6 1.4 1245 122
## [4,] 11.3 13.2 3.9 742.2 27.8 -5.6 33.4 18.9 1.4 20.6 1.4 1245 122
## [5,] 13.1 13.8 3.9 767.0 30.3 -4.3 34.6 21.0 4.5 22.8 2.9 1183 113
## [6,] 13.1 13.8 3.9 767.0 30.3 -4.3 34.6 21.0 4.5 22.8 2.9 1183 113
## pdry pseas pwetq pdryq pwarmq pcoldq
## [1,] 79 9 304 252 304 259
## [2,] 83 11 345 264 329 264
## [3,] 83 11 345 264 329 264
## [4,] 83 11 345 264 329 264
## [5,] 84 10 328 257 318 264
## [6,] 84 10 328 257 318 264nrow(fire_climate)## [1] 1188firenc <- cbind(firenc, fire_climate)
head(firenc)## Simple feature collection with 6 features and 36 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -81.108 ymin: 36.421 xmax: -79.374 ymax: 36.543
## Geodetic CRS: NAD83
## AREA PERIMETER FIRE_ FIRE_ID LAT LONG DATE JULIAN GMT TEMP
## 49833 0 0 49833 259985 36.543 -79.374 2010-08-12 224 1822 322.1
## 50331 0 0 50331 337411 36.490 -81.108 2010-10-20 293 1840 315.8
## 50332 0 0 50332 337410 36.490 -81.108 2010-10-20 293 1840 315.8
## 50350 0 0 50350 337409 36.489 -81.107 2010-10-20 293 1840 315.8
## 51019 0 0 51019 470917 36.421 -80.686 2010-11-22 326 1603 324.2
## 51020 0 0 51020 470918 36.421 -80.686 2010-11-22 326 1604 324.2
## SPIX TPIX SRC SAT_SRC CONF FRP ID mat mdr iso tseas
## 49833 1.0 1 gsfc_drl A 51 12.2 north carolina 14.1 13.5 3.8 787.9
## 50331 1.1 1 gsfc_drl A 64 21.9 north carolina 11.3 13.2 3.9 742.2
## 50332 1.1 1 ssec A 64 21.9 north carolina 11.3 13.2 3.9 742.2
## 50350 1.1 1 gsfc A 64 21.9 north carolina 11.3 13.2 3.9 742.2
## 51019 1.1 1 ssec T 77 31.7 north carolina 13.1 13.8 3.9 767.0
## 51020 1.1 1 gsfc_drl T 77 31.7 north carolina 13.1 13.8 3.9 767.0
## tmax tmin tar twetq tdryq twarmq tcoldq ap pwet pdry pseas pwetq pdryq
## 49833 31.7 -3.3 35.0 24.1 5.4 24.1 3.7 1112 111 79 9 304 252
## 50331 27.8 -5.6 33.4 18.9 1.4 20.6 1.4 1245 122 83 11 345 264
## 50332 27.8 -5.6 33.4 18.9 1.4 20.6 1.4 1245 122 83 11 345 264
## 50350 27.8 -5.6 33.4 18.9 1.4 20.6 1.4 1245 122 83 11 345 264
## 51019 30.3 -4.3 34.6 21.0 4.5 22.8 2.9 1183 113 84 10 328 257
## 51020 30.3 -4.3 34.6 21.0 4.5 22.8 2.9 1183 113 84 10 328 257
## pwarmq pcoldq geometry
## 49833 304 259 POINT (-79.374 36.543)
## 50331 329 264 POINT (-81.108 36.49)
## 50332 329 264 POINT (-81.108 36.49)
## 50350 329 264 POINT (-81.107 36.489)
## 51019 318 264 POINT (-80.686 36.421)
## 51020 318 264 POINT (-80.686 36.421)plot(TEMP ~ ap, data=firenc) 
Here is a quick demo of a new R package called
leaflet
https://rstudio.github.io/leaflet
library(leaflet)
mapStates = map("state", fill = TRUE, plot = FALSE)
leaflet(data = mapStates) %>% addTiles() %>%
addPolygons(fillColor = topo.colors(10, alpha = NULL), stroke = FALSE)m <- leaflet(data = firenc) %>% addTiles() %>%
addCircleMarkers(radius = 2, label = ~as.character(firenc$TEMP))## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'm %>% addProviderTiles(providers$Esri.NatGeoWorldMap)m = leaflet() %>% addTiles()
df = data.frame(
lat = rnorm(100),
lng = rnorm(100),
size = runif(100, 5, 20),
color = sample(colors(), 100)
)
m = leaflet(df) %>% addTiles()
m %>% addCircleMarkers(radius = ~size, color = ~color, fill = FALSE)## Assuming "lng" and "lat" are longitude and latitude, respectivelym %>% addCircleMarkers(radius = runif(100, 4, 10), color = c('red'))## Assuming "lng" and "lat" are longitude and latitude, respectivelym %>% addProviderTiles(providers$Esri.NatGeoWorldMap)