# Distance to a lake How far are you from a lake?

Living in Wisconsin, I wanted to know. 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)
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))