Living 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"),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')