# Create density polygons from spatial points in R

I have a dataframe with a bunch of spatial points representing flooding occurrences. I want to know which areas have the highest density of points so we can prioritize areas for management. In ArcGIS Pro you can do this by changing the symbology to create a heatmap.

It looks like the `ggmap` library has a `stat_density_2d()` function which might do the trick. The thing is I don't just want a map/plot, I want to export the polygons/contours as a shapefile to use in a series of analyses.

What I have so far:

``````#loading libraries
library(ggplot2)
library(ggmap)
library(sf)
library(dplyr)

#creating polygon
polygon <- list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1),ncol=2, byrow=T))
polygon <- sf::st_polygon(polygon) # Create an sf polygon

# Sample 50 random points within the polygon
points <- st_sample(polygon, size=50)
st_crs(points) <- "+proj=longlat +ellps=WGS84 +datum=WGS84" #change to lon/lat
points <- as.data.frame(points) #convert to dataframe
points\$lon <- unlist(map(points\$geometry,1)) #adding lon column
points\$lat <- unlist(map(points\$geometry,2)) #adding lat column

#Plot using geom_density_2d()
ggplot(points, aes(x = lon, y = lat)) +
geom_point() +
geom_density_2d()

``````

The issue is I want those contours as polygons that I can export and use in other analyses

• Note contours aren't always polygons on a map, since they can run off the edge. Nov 9, 2022 at 16:34

Starting with random points as a spatial data frame:

``````points <- st_sample(polygon, size=50)
st_crs(points) <- "+proj=longlat +ellps=WGS84 +datum=WGS84" #change to lon/lat
``````

get the X-Y coords (lon-lat):

``````> xy = st_coordinates(points)
X        Y
1 1.285699 1.716926
2 1.057654 1.937922
3 1.603304 1.753047
``````

Use `MASS::kde2d` which is what `geom_density_2d` uses, it says in the docs:

``````> de = MASS::kde2d(xy[,1], xy[,2])
> image(de)
``````

That can be converted to a raster for export to other systems. If you want actual contour lines, these are going to be dependent on R's contouring algorithm. You can get them for example with:

``````cl = contourLines(de)
``````

This returns a list of x-y lists which you can plot on top of your raster:

``````lapply(cl, lines)
``````

The value of each line is stored in the level element alongside the X and Y for that line:

``````> cl[[1]]\$level
[1] 0.2
``````

so you could write a nice contour plot function that labelled the lines. This probably exists somewhere. You could also convert this to a spatial lines object with the level in the attribute table for export, but since you've not been precise about where you want to export this to I'll leave it for now.

Edit: here's a one-liner (split for readability) that converts the `cl` list of x,y, value to a spatial lines data frame with value:

``````cllines = do.call(rbind,
Map(function(x){
st_as_sf(
data.frame(
Z=x\$level,
geometry=st_sfc(
st_linestring(cbind(x\$x, x\$y))
)
)
)},cl
))

``````

`plot(cllines)` will then give you a coloured lines plot with scale bar.

I forgot the BIG WARNING. Don't do this with real lat-long data, because now your kernel smoothing isn't the same distance in latitude and longitude. Project your data to a locally cartesian system and work with that.

• On second thought I feel I have misunderstood the meaning of "areas" in original question; deleting my answer, as yours is clearly better. respect. Nov 9, 2022 at 21:11
• Thanks so much for the great answer! To convert `cl` to a spatial lines object. Would I need to write a function that converts each item in the list? Nov 15, 2022 at 18:40
• Have now added a function to convert the lines list to spatial lines data frame with height value... it uses `Map` to apply a conversion function to each element to make a data frame of one row and then rbind them all together.. Nov 15, 2022 at 21:29

You can extract the underlying data from your ggplot2 after transformation and convert it to a sf polygon. There is a slight gotcha, you would need to use `geom_density2d_filled()` to get the info for building polygons:

``````#loading libraries
library(ggplot2)
library(sf)
library(dplyr)
library(purrr)

#creating polygon
polygon <- list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1),ncol=2, byrow=T))
polygon <- sf::st_polygon(polygon) # Create an sf polygon

# Sample 50 random points within the polygon
set.seed(1234)

points <- st_sample(polygon, size=50)
st_crs(points) <- "+proj=longlat +ellps=WGS84 +datum=WGS84" #change to lon/lat

points_sf <- points

points <- as.data.frame(points) #convert to dataframe
points\$lon <- unlist(map(points\$geometry,1)) #adding lon column
points\$lat <- unlist(map(points\$geometry,2)) #adding lat column

#Plot using geom_density_2d_filled()
p <- ggplot(points, aes(x = lon, y = lat)) +
# This would create filled polygons that we use for create our polygon
geom_density2d_filled() +
geom_point() +
geom_density_2d()

p
``````

Now we extract the information of the filled contour layer (the first layer of our plot) and we covert it to sf polygons:

``````# Extract underlying data from ggplot2, first layer
data2d <- layer_data(p, i = 1)

#>        fill      level        x        y piece  group subgroup level_low
#> 1 #440154FF (0.0, 0.2] 1.992150 1.926400     1 -1-001        1         0
#> 2 #440154FF (0.0, 0.2] 1.982225 1.926400     1 -1-001        1         0
#> 3 #440154FF (0.0, 0.2] 1.972299 1.926400     1 -1-001        1         0
#> 4 #440154FF (0.0, 0.2] 1.969325 1.926400     1 -1-001        1         0
#> 5 #440154FF (0.0, 0.2] 1.972299 1.917381     1 -1-001        1         0
#> 6 #440154FF (0.0, 0.2] 1.972367 1.917182     1 -1-001        1         0
#>   level_high level_mid nlevel  n PANEL colour linewidth linetype alpha
#> 1        0.2       0.1    0.1 50     1     NA       0.5        1    NA
#> 2        0.2       0.1    0.1 50     1     NA       0.5        1    NA
#> 3        0.2       0.1    0.1 50     1     NA       0.5        1    NA
#> 4        0.2       0.1    0.1 50     1     NA       0.5        1    NA
#> 5        0.2       0.1    0.1 50     1     NA       0.5        1    NA
#> 6        0.2       0.1    0.1 50     1     NA       0.5        1    NA

# Id of polygon
data2d\$pol <- paste0(data2d\$group, "_", data2d\$subgroup)
ids <- unique(data2d\$pol)

# Split and create polygons based on the id
pols <- lapply(ids, function(x){
topol <- data2d[data2d\$pol == x, ]

closepol <- rbind(topol, topol[1, ])

pol <- st_polygon(list(as.matrix(closepol[,c("x", "y")])))

df <- unique(topol[, grepl("level", names(topol))])

tofeatures <- st_as_sf(df, geometry=st_sfc(pol))

return(tofeatures)
})

final_pols <- do.call(rbind, pols)

# And force a crs, since we lost that on the process
st_crs(final_pols) <- st_crs(points_sf)

plot(final_pols["level"], axes = TRUE)
``````

``````# Check

ggplot(final_pols) +
geom_sf(aes(fill=level)) +
geom_sf(data = points_sf)
``````

Created on 2022-11-10 with reprex v2.0.2

• Note that here I am using ggplot2 for processing data, that it is not their (primary) intended application, but still hope this could help. Cheers Nov 10, 2022 at 18:21
• Thanks for your help this was really great! Dec 12, 2022 at 16:24