2

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

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

2 Answers 2

3

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)
> head(xy)
         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)

enter image description here

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.

enter image description here

3
  • 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?
    – seak23
    Nov 15, 2022 at 18:40
  • 1
    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..
    – Spacedman
    Nov 15, 2022 at 21:29
2

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) 

head(data2d)
#>        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")])))
  
  # Add features
  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

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
    – dieghernan
    Nov 10, 2022 at 18:21
  • Thanks for your help this was really great!
    – seak23
    Dec 12, 2022 at 16:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.