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
rnorm to generate 20000 coordinate pairs randomly. You could also use
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
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.