2

I Have a shapefile of polygons and another one of points that are distributed over the polygons. I would like to create a kernel density estimate for each polygon based on the points it contains. unfortunately I was only able to create squared KDEs with the kde2d function from the MASS package. I would like the KDEs to be shaped as the polygons. Any suggestions?

 kde1 <- kde2d(poly$X, poly$Y, n=100,)

enter image description here

Ege Rubak
  • 4,347
  • 1
  • 10
  • 18
Naama
  • 39
  • 4

2 Answers2

3

You can use the spatstat package for this. Here is an example of reading in a shapefile from sf, generating random points and run kernel density estimation of the intensity of points (points per unit area):

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
nc <- st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source `/usr/lib/R/site-library/sf/shape/nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> geographic CRS: NAD27
nc_flat <- st_transform(nc, crs = 26917)
W <- as.owin(nc_flat$geometry[1]) # First county of North Carolina data set in spatstat format

library(spatstat)
X <- runifpoint(100, win = W)
plot(X, "Random points")


D <- density(X)
plot(D, main = "KDE")

Ege Rubak
  • 4,347
  • 1
  • 10
  • 18
  • Great, thank you! But I failed to make it work with my own points shapefile (and not the random points as you suggested). I keep getting the " argument 'x' must be numeric" error, and I'm not sure what should I do to make the shapefile suitable for the "density" function. Do you have any suggestions? – Naama Feb 01 '21 at 15:37
  • OK! I managed to use my own points by using the 'ppp' function from the spatstat package. ` C <- as.owin(polygon$geometry[p]) p<- ppp(points$X,points$Y, window = C) D <- density(p)` – Naama Feb 01 '21 at 16:59
0

OK! I managed to use my own points by using the 'ppp' function from the spatstat package.

  C <- as.owin(polygon$geometry[n])

  p<- ppp(points$X,points$Y, window = C)

  D <- density(p)
[enter image description here][1]


  [1]: https://i.stack.imgur.com/YZN0V.png
Naama
  • 39
  • 4
  • The command `D <- density(p)` is using the function `density.ppp` which is the method for the generic `density` for point patterns (class `ppp`). See the help for `density.ppp` to find out how to control the amount of smoothing, etc. – Adrian Baddeley Feb 03 '21 at 08:04
  • @AdrianBaddeley, Thank you. For some reason, when I use both the `density` and the `density.ppp`, the output is different than the one described in the function's document. Instead of a list containing 'x', 'y', 'bw' and 'n', i get a list containing 'v', 'xrange', 'yrange' and a few more, but no kernel values i can use. any idea why this happens? – Naama Feb 03 '21 at 18:46