Distance to a lake

wilakedistance2Living in Wisconsin, I was curious how far from a lake I was at any given point in the state. There’s a lot of lakes in Wisconsin, but it was hard to visualize their distribution.

Step 1) Download the high-resolution National Hydrography Dataset for Wisconsin (I host this locally as a shapefile). I only kept lakes with a surface area larger than 4 ha, and only polygons designated as lakes ‘Lake/pond; a standing body of water with a predominantly natural shoreline surrounded by land.’

nhd = readOGR(dsn = instertlocaldirectory,layer = 'WI')
areaX = nhd@data$AreaSqKm > 0.04 # Lakes > 4 ha
lakeX = nhd@data$FType == 390 # Only lakes
nhd2 = nhd[areaX & lakeX,]

The NHD data is presented in lat/long. Since I want to measure distances, I needed to project it. Wisconsin is in UTM Zone 15.

nhdUTM = spTransform(nhd2,CRSobj = "+proj=utm +zone=15 ellps=WGS84")

Step 2) Rasterize lake data onto a grid. The rasterize function (package:: raster) transfer values associated with ‘object’ type spatial data (points, lines, polygons) to raster cells.

emptyRaster = raster(ncol=5000, nrow=4280) #approximately 100 m grid for WI 
emptyRaster = extent(nhdUTM) #add spatial extent
lakeRaster = rasterize(nhdUTM,emptyRaster,field = 1, background = 0)
lakeRaster_class = reclassify(lakeRaster,cbind(0,NA))

Step 3) Compute distance from each grid cell to closest cell with value “1”. The raster::distance function does the heavy lifting here.

lakeDistance = distance(lakeRaster_class)

Step 4) Crop our data to just the state of Wisconsin. I can get a state map from the maps package in R.

library(maps); library(mapStats) 
usMap <- readOGR(system.file("shapes/usMap.shp", package="mapStats")[1],layer='usMap')
stMap = usMap[usMap@data$STATE_ABBR==state,]
stMap = spTransform(stMap,crs(nhdUTM))

Crop lake distance grid with state map:

distanceWI = rasterize(stMap,lakeDistance,mask=T)

Step 5) Plot the results! (in miles for convenience)

par(mar=c(2,0,2,0),mgp=c(1.5,0.5,0))
plot(distanceWI*0.62/1000,breaks = c(0,2,5,10,20),axes=FALSE, box=FALSE,
 col = colorRampPalette(c('azure2', "navy"))( 5),
 main = 'Distance from a lake in Wisconsin',
 legend.args=list(text='Distance (miles)', side=4, font=2, line=2, cex=0.8),
 axis.args=list(at=c(0,2,5,10,20), labels=c(0,2,5,10,'>20'), cex.axis=1))
plot(stMap,add=T,border='grey30')
Advertisements

Extracting polygon areas in R

Suppose I have a shapefile (lakes.shp) which consists of 5 polygons (Figure 1a).

library(raster)
lakes = shapefile('lakes.shp')
plot(lakes, col='dodgerblue',main='a) lakes',lty=0)
Figure 1: a) Lake polygons, b) Lake polygons with 30 m buffer
Figure 1: a) Lake polygons, b) Lake polygons with 30 m buffer

The shapefile is stored as a SpatialPolygonsDataFrame. The data underlying these polygons can be viewed with:

lakes@data

There is a column called SHAPE_Area, which is the area of each polygon. Now let’s say, we want to create a buffer zone around the lake (Figure 1b).

lakebuffer = buffer(lakes, width=30, dissolve=F) 
plot(buffer_ring,col='darkred',main = 'b) 30 m buffer',lty=0)

Great! We can still see the underlying data with

buffer_ring@data

However, the SHAPE_Area shown are the original lake areas. To find the true polygon area, we must dig into the polygon layer: buffer_ring@polygons. But each polygon is an element of a list, and the slot with “area” is deeply nested. Instead of looping through each polygon, here’s an easy way to extract polygon areas:

bufAreas = sapply(slot(lakebuffer, "polygons"), slot, "area")

For the area of only the buffer:

lakeAreas = sapply(slot(lakes, "polygons"), slot, "area")
output = data.frame(bufferArea = bufAreas - lakeAreas)