# Kernel Density Estimation in R

For a recent project I needed to run a kernel density estimation in R, turning GPS points into a raster of point densities. Below is how I accomplished that.

## Load the needed libraries

`library(sf)`

library(raster)

library(ggspatial)

If you don't have these packages you need to run `install.packages(c("sf", "raster"))`

before they can be loaded.

## Create XY Points

For the purposes of this exercise, and because the original dataset I worked with is not available to the public, we will not be loading a shapefile but instead making one out of randomly generated X and Y coordinates. We will create UTM coordinates in Zone 13. You could also create Latitude and Longitude coordinates if you wanted. Here are the bounds of the coordinates:

**Northing max:** 4232607

**Northing min:** 4220932

**Easting max:** 187973

**Easting min:** 173118

We'll use `rnorm`

to generate 20000 coordinate pairs randomly. You could also use `runif`

but `rnorm`

creates areas with higher densities so our plot is not uniform.

`coord<-tibble(`

northing = rnorm(10000, 4220105, 3589.521),

easting = rnorm(10000, 185310.8, 12709.06)

)

Now we can plot the results to verify the points were generated correctly.

`coord%>%`

ggplot(aes(northing, easting))+

geom_point(color = "steelblue", size = 3, alpha = 0.35)+

theme_minimal()

Everything looks good. Now we need to convert the xy coordinates to a shapefile.

## Converting xy coordinates to a shapefile

Now that we have xy pairs we can easily convert the xy pairs to a shapefile and give the shapefile a projection. If you need more info on projections you can find it here

`sf_obj<- st_as_sf(coord, coords = c("easting", "northing"), crs = "+proj=utm +zone=13 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")`

sf_obj

Output:

`Simple feature collection with 5000 features and 0 fields`

geometry type: POINT

dimension: XY

bbox: xmin: 173118.4 ymin: 4220933 xmax: 187971.9 ymax: 4232603

epsg (SRID): 26913

proj4string: +proj=utm +zone=13 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

# A tibble: 5,000 x 1

geometry

<POINT [m]>

1 (176681.1 4230659)

2 (174927.8 4222673)

3 (178167.1 4229261)

4 (186605.6 4230476)

5 (176185.4 4227015)

6 (184934.5 4222134)

7 (186177.4 4231785)

8 (183417 4231574)

9 (187055.3 4222734)

10 (174989 4230511)

# ... with 4,990 more rows

Let's take a look at our new shapefile.

`ggplot()+`

geom_sf(data = sf_obj, color = "steelblue", size = 3, alpha = 0.3)+

theme_minimal()

The plot looks very similar to the xy plot (they should be the same), but now the points are projected. If you look at the lines behind the points you see that they are slightly askew.

## Make a empty raster

Now we will create an empty raster that we will use to calculate the kernel density from. The resolution we set will be the resolution of the kernel density. We will use the `raster()`

function from the raster package to make the raster.

`empty_kernel_grid<-raster(ext = extent(sf_obj), resolution = 200, crs = st_crs(sf_obj))`

Next we will take our points that we created and use the `fun = 'count'`

parameter in the `rasterize`

function from the to count all the points that fall within a given raster cell. We will also tell rasterize to use the `empty_kernel_grid`

raster for the bounds of the raster.

`kernel_density<-rasterize(coordinates(as_Spatial(sf_obj)), empty_kernel_grid, fun='count', background = 0)`

To plot a raster you with ggplot, you first must convert the raster to points with rasterToPoints.

`library(RColorBrewer)`

point_raster<-as_tibble(rasterToPoints(kernel_density))

ggplot(point_raster, aes(x = x, y=y))+

geom_raster(aes(fill = layer))+

theme_minimal()+

scale_fill_distiller(palette = "Spectral")

Even though we had a mass to start, now we see that certain 200m areas are mor dense than other.