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')