1.1 Calculate density (km/km2) of linear features

This is the first somewhat tricky task. We’ll do this for a sample portion of the map first to speed up processing time. st_crop crops the shapefile to the listed extent:

require(raster)
BEAD.crop <- st_crop(BEAD.lines, extent(-1253023, -1227564, 2546649, 2578169))
plot(BEAD.crop[,"Class"])

BEAD.crop.sp <- as(BEAD.crop, "Spatial")

as(,"Spatial") converts an sf object to a spatial object. Alternatively, you can use as.Spatial().

We use a multi-step process to calculate the density of linear features:

1. Make an empty raster with the same extent as the linear features shapefile
2. Make a shapefile of the empty raster
3. Cut the linear features at the transitions between pixels using the new shapefile
4. Populate the empty raster with the sum of the lengths of the linear features that intersect it.

Steps 2 and 3 are important for the accuracy of the density calculation when using this method. The method used here measures the length of all linear features that intersect with a pixel and adds them together. That means if a road intersects with a pixel, it will measure the length of the entire road, not just the part within that pixel. Steps 2 and 3 cut the linear features at the borders between pixels, so that when you measure the lengths in step 4 you count only the portion of each linear feature that is inside that particular pixel.

require(plyr)
require(rgeos)

First we make an empty raster. It should have the same extent and projection as the cropped shapefile:

empty.raster <- raster(extent(BEAD.crop), resolution = 3000, crs = projection(BEAD.crop))

Next, make a shapefile from the empty raster that looks like a grid, i.e. the outline of each pixel. rasterToPolygons(raster package) makes a shapefile of the outline of each pixel.

density.poly <- rasterToPolygons(empty.raster)
plot(density.poly)

The tricky step is to “intersect” the lines with his grid, which cuts the linear features at the edges.gIntersection(rgeos package) is a function for determining the intersection between two given geometries (in this case, density.poly and BEAD.crop.sp). Note, the byid = TRUE in the gIntersection function makes sure that the clipping occurs for EACH polygon. The later stuff converts it back to a simple feature.

require(magrittr)
lines.clipped <- gIntersection(density.poly, BEAD.crop.sp, byid = TRUE) %>% 
  st_as_sf %>% mutate(length = st_length(geometry))

We now have a simple feature with segments cut into the new raster with a computed “length” column:

lines.clipped

## Simple feature collection with 173 features and 1 field
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -1253023 ymin: 2546649 xmax: -1229023 ymax: 2578169
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## First 10 features:
##                          geometry       length
## 1  LINESTRING (-1252069 257794...  235.85716 m
## 2  LINESTRING (-1253023 257602...  216.56771 m
## 3  LINESTRING (-1242183 257537...  222.31043 m
## 4  LINESTRING (-1243536 257816... 3130.52798 m
## 5  LINESTRING (-1243745 257816...    7.74042 m
## 6  LINESTRING (-1241023 257809... 2954.10981 m
## 7  LINESTRING (-1242183 257537...  219.81122 m
## 8  LINESTRING (-1240989 257816...   85.83307 m
## 9  LINESTRING (-1238037 257556...  412.51364 m
## 10 LINESTRING (-1235023 257568...   97.96807 m

plot(density.poly, border = "grey")
plot(lines.clipped, add = TRUE, lwd = 2)

Finally, we populate the empty raster with the sum of the lengths of the linear features that intersect it using the rasterize command. The.

lines.clipped$length <- as.numeric(lines.clipped$length)/1e3
density.raster <- rasterize(as(lines.clipped, "Spatial"),
                                   empty.raster, 
                                   field = "length", fun="sum")/9

The rasterize is pretty slow Extra credit: Is there a way to do this very quickly with fasterize? I’m not sure!

Note - the first line converts the distance to km. The ./9 in the second line takes into account the fact that our raster is 3km x 3km, so to get to km/km^2 we divide by 9.

Plot the linear features over your new density raster

palette(topo.colors(100))
plot(density.raster/9)
plot(BEAD.crop[,"Class"], add = TRUE, lwd = 2, pal = rainbow)

OR, you can use mapview for a nicer visualization against a map:

m1 <- mapview(density.raster,  map.types = "Esri.WorldImagery")
mapview(BEAD.crop[,"Class"], map = m1@map)

The mapview features are great for confirming that these spatial analyses are working, as you can zoom in to any particular pixel, and scroll over any particular line segment.

Some key functions used in this code:

• rasterToPolygons - in raster
• rasterize - in raster
• st_length - in sf
• Several times we had to convert between “Spatial” and “simple feature” type objects using st_as_sf (converts spatial to sf) and as(., "Spatial") - (sf to spatial).

getLineDensity <- function(lines.sf, resolution, plotme = TRUE){
  pcks <- c("plyr","rgeos","sp","sf","raster")
  a <- sapply(pcks, require, character = TRUE)
  
  empty.raster <- raster(extent(lines.sf), 
                         resolution = resolution, 
                         crs = projection(lines.sf))
  
  density.poly <- rasterToPolygons(empty.raster)
  
  lines.sp <- as(lines.sf, "Spatial")
    
  lines.clipped <- gIntersection(density.poly, lines.sp, byid = TRUE) %>% 
    st_as_sf %>% mutate(length = st_length(geometry))
  
  lines.clipped$length <- as.numeric(lines.clipped$length)/1e3
  
  density.raster <- rasterize(as(lines.clipped, "Spatial"),
                              empty.raster, field = "length",
                              fun="sum")/((resolution/1e3)^2)
  if(plotme){
    plot(density.raster)
    plot(st_geometry(lines.sf), add = TRUE, lwd = 2, pal = rainbow)
  }
  return(density.raster)
}

line.density1km <- getLineDensity(BEAD.crop, 1e3)

system.time(line.density10km <- getLineDensity(subset(BEAD.lines, Class == "Road"), 10e3))

##    user  system elapsed 
##   17.31    0.06   17.89

system.time(line.density10km <- getLineDensity(BEAD.lines, 10e3))

##    user  system elapsed 
##   95.14    0.39   98.79

m1 <- mapview(line.density10km,  map.types = "Esri.WorldImagery")
mapview(BEAD.lines[,"Class"], map = m1@map)

1.2 Distance to linear elements

To obtain a raster of distances from nearest linear elements, we:

1. find the centroids of the desired raster with st_centroid
2. find the nearest distances between that set of centroids and the linear elements with st_distance

Find centroids of the raster:

empty.raster <- raster(extent(BEAD.crop), resolution = 3000, crs = projection(BEAD.crop))
distance.poly <- rasterToPolygons(empty.raster) 
distance.poly.sf <- st_as_sf(distance.poly)

Now we can compute centroids:

distance.centroid <- st_centroid(distance.poly.sf)
plot(distance.poly)
plot(st_geometry(distance.centroid), add = TRUE)

There is a slight hiccup before the next step: for some reason the conversion from spatial slightly changes the projection definition, which causes problems later …

st_crs(distance.centroid)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

st_crs(BEAD.crop)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

distance.centroid <- st_transform(distance.centroid, st_crs(BEAD.crop))
st_crs(distance.centroid) == st_crs(BEAD.crop)

## [1] TRUE

Obtain distances:

distances <- st_distance(distance.centroid, BEAD.crop)

This generates a large matrix of all pairwise distances:

nrow(distance.centroid)

## [1] 88

nrow(BEAD.crop)

## [1] 116

dim(distances)

## [1]  88 116

We want the closest distance, which we obtain as:

distance.min <- apply(distances, 1, min)

We can now populate the empty raster with these values:

distance.raster <- setValues(empty.raster, distance.min)

Let’s check this:

plot(distance.raster)
plot(st_geometry(BEAD.crop), add = TRUE)

OR - with mapview:

m1 <- mapview(distance.raster, map.types = "Esri.WorldImagery")
mapview(BEAD.crop, m1@map)

getDistanceToLines <- function(lines.sf, resolution, plotme = TRUE){

  pcks <- c("plyr","rgeos","sp","sf","raster", "fasterize")
  a <- sapply(pcks, require, character = TRUE)
  
  xy <- st_coordinates(lines.sf)[,1:2]
  mcp <- xy[chull(xy),] %>% rbind(xy[chull(xy)[1],]) %>% list %>% st_polygon %>% st_sfc %>% st_sf
  empty.raster <- fasterize(mcp, raster(extent(lines.sf), resolution = resolution, crs = projection(lines.sf)))

  # cut out portions of raster outside MCP
  distance.centroid <- suppressWarnings(rasterToPolygons(empty.raster) %>% st_as_sf %>% st_centroid %>% st_transform(st_crs(lines.sf)))
  
  distance.matrix <- st_distance(distance.centroid, lines.sf)
  distance.min <- apply(distance.matrix, 1, min)
  
  raster.vals <- getValues(empty.raster)
  raster.vals[!is.na(raster.vals)] <- distance.min

  distance.raster <- setValues(empty.raster, raster.vals)
  
  if(plotme){
    plot(distance.raster)
    plot(st_geometry(lines.sf), add = TRUE)
  }
  
  return(distance.raster)
}

system.time(d.to.line_1km <- getDistanceToLines(BEAD.crop, 1e3))

##    user  system elapsed 
##    0.66    0.00    0.70

system.time(d.to.road_10km <- getDistanceToLines(subset(BEAD.lines, Class == "Road"), 10e3))

##    user  system elapsed 
##    0.55    0.00    0.58

system.time(d.to.lines_10km <- getDistanceToLines(BEAD.lines, 5e3))

##    user  system elapsed 
##   15.76    0.13   16.31

2.1 Calculate density of polygonal features

We will take the polygon data and compute the density into a 1km raster. In the code below, there are 5 steps:

1. Create an empty 1x1 km raster:

empty.raster <- raster(extent(BEAD.polys.crop), resolution = 1e3, crs = projection(BEAD.polys.crop))

2. Convert to a polygon data type. Note that rasterToPolygons makes a Spatial data object, NOT a simple feature… this is another function that some time soon (if not already) will be doable with sf:

polys.grid.st <- rasterToPolygons(empty.raster) 

3. Obtain the interection between the two poygon data sets. Note, we have to convert the BEAD polygon to Spatial as well. Note also, that the byid = TRUE is necessary to get the very specific slicing of the two polygon data sets:

require(rgeos)
polys.raster.intersection <- gIntersection(polys.grid.st, as(BEAD.polys.crop, "Spatial"), byid = TRUE)
plot(polys.raster.intersection)

4. Importantlty, we need to extract the index of the appropriate raster polygon, to correctly populate the raster. Note that the names of the ensuing intersected polygons combine the indices of each of the polygons in the original dataset:

head(names(polys.raster.intersection))

## [1] "17 60"  "17 182" "18 182" "19 182" "34 185" "35 61"

We need the first number in each pair. The following bit of code is one way to do that:

index <- sapply(strsplit(names(polys.raster.intersection), " "), head, 1) %>% as.numeric 

5. Extract the area of each of the polygons. I find this is easier switching back to simple features:

areas <- st_area(st_as_sf(polys.raster.intersection))

6. Populate the raster - but carefully, because several of the polygons have multiple area values to fill. So first, we need to compute the total area in each raster block:

area.final <- tapply(areas, index, sum)
index.final <- unique(index)

Now we know exactly where to place these.

empty.raster[index.final] <- area.final/1e6
plot(empty.raster)
plot(BEAD.polys.crop, add = TRUE)

2.2 In a single function

All together as a function. Note that I added an “overlapping” tag. If the polygons are NOT overlapping, the function is much faster because it can “collapse” all of the polygons into a single object with the gUnaryUnion function:

getPolygonDensity <- function(polys.sf, resolution, plotme = TRUE, overlapping = FALSE){
  pcks <- c("plyr","rgeos","sp","sf","raster","fasterize")
  
  a <- sapply(pcks, require, character = TRUE)
  density.raster <- raster(extent(polys.sf), 
                         resolution = resolution, 
                         crs = projection(polys.sf))
  
  raster.poly <- rasterToPolygons(density.raster)
  polys.sp <- as(polys.sf, "Spatial") 
  if(!overlapping) polys.sp <- gUnaryUnion(polys.sp)
  
  polys.raster.intersection <- gIntersection(raster.poly, polys.sp, byid = TRUE)
  
  pri.sf <- st_as_sf(polys.raster.intersection)  
  pri.sf$area <- st_area(pri.sf)
  
  index <- sapply(strsplit(names(polys.raster.intersection), " "), head, 1) %>% as.numeric 
  areas <- st_area(st_as_sf(polys.raster.intersection))
  area.final <- tapply(areas, index, sum)
  index.final <- unique(index)
  density.raster[index.final] <- area.final/1e6

  if(plotme){
    plot(density.raster)
    plot(st_geometry(polys.sf), add = TRUE)
  }
  return(density.raster)
}

Some examples - at 1km:

system.time(getPolygonDensity(BEAD.polys.crop, 1000))

##    user  system elapsed 
##    0.40    0.00    0.45

At 0.5 km:

system.time(getPolygonDensity(BEAD.polys.crop, 500))

##    user  system elapsed 
##    0.98    0.02    1.00

At 5 km for entire range, comparing overlapping and non-overlapping:

system.time(a <- getPolygonDensity(BEAD.polys, 5000, overlapping = TRUE))

##    user  system elapsed 
##    8.66    0.05    8.92

system.time(b <- getPolygonDensity(BEAD.polys, 5000, overlapping = FALSE))

##    user  system elapsed 
##    5.02    0.04    5.15