R - 如何在特定轮廓内找到点

I am creating density plots with kde2d (MASS) on lat and lon data. I would like to know which points from the original data are within a specific contour.

我正在使用lat和lon数据上的kde2d（MASS）创建密度图。我想知道原始数据中哪些点在特定轮廓内。

I create 90% and 50% contours using two approaches. I want to know which points are within the 90% contour and which points are within the 50% contour. The points in the 90% contour will contain all of those within the 50% contour. The final step is to find the points within the 90% contour that are not within the 50% contour (I do not necessarily need help with this step).

我使用两种方法创建90％和50％的轮廓。我想知道哪些点在90％轮廓内，哪些点在50％轮廓内。 90％轮廓中的点将包含50％轮廓内的所有点。最后一步是找到90％轮廓内不在50％轮廓内的点（我不一定需要这个步骤的帮助）。

# bw = data of 2 cols (lat and lon) and 363 rows
# two versions to do this: 
# would ideally like to use the second version (with ggplot2)

# version 1 (without ggplot2) 
library(MASS)
x <- bw$lon
y <- bw$lat
dens <- kde2d(x, y, n=200)

# the contours to plot
prob <- c(0.9, 0.5)
dx <- diff(dens$x[1:2])
dy <- diff(dens$y[1:2])
sz <- sort(dens$z)
c1 <- cumsum(sz) * dx * dy 
levels <- sapply(prob, function(x) { 
    approx(c1, sz, xout = 1 - x)$y
})
plot(x,y)
contour(dens, levels=levels, labels=prob, add=T)

And here is version 2 - using ggplot2. I would ideally like to use this version to find the points within the 90% and 50% contours.

这是版本2 - 使用ggplot2。我最好使用这个版本来找到90％和50％轮廓内的点。

# version 2 (with ggplot2)
getLevel <- function(x,y,prob) { 
    kk <- MASS::kde2d(x,y)
    dx <- diff(kk$x[1:2])
    dy <- diff(kk$y[1:2])
    sz <- sort(kk$z)
    c1 <- cumsum(sz) * dx * dy
    approx(c1, sz, xout = 1 - prob)$y
}

# 90 and 50% contours
L90 <- getLevel(bw$lon, bw$lat, 0.9)
L50 <- getLevel(bw$lon, bw$lat, 0.5)

kk <- MASS::kde2d(bw$lon, bw$lat)
dimnames(kk$z) <- list(kk$x, kk$y)
dc <- melt(kk$z)

p <- ggplot(dc, aes(x=Var1, y=Var2)) + geom_tile(aes(fill=value)) 
+ geom_contour(aes(z=value), breaks=L90, colour="red")
+ geom_contour(aes(z=value), breaks=L50, color="yellow")
+ ggtitle("90 (red) and 50 (yellow) contours of BW")

I create the plots with all of the lat and lon points plotted and 90% and 50% contours. I simply want to know how to extract the exact points that are within the 90% and 50% contours.

我创建了绘制了所有lat和lon点以及90％和50％轮廓的图。我只是想知道如何提取90％和50％轮廓内的确切点。

I have tried to find the z values (the elevation of the density plots from kde2d) that are associated with each row of lat and lon values but had no luck. I was also thinking I could add an ID column to the data to label each row and then somehow transfer that over after using melt(). Then I could simply subset the data that has values of z that match each contour I want and see which lat and lon they are compared to the original BW data based on the ID column.

我试图找到与每一行lat和lon值相关联但没有运气的z值（来自kde2d的密度图的高程）。我还在想我可以在数据中添加一个ID列来标记每一行，然后在使用melt（）之后以某种方式将其转移。然后，我可以简单地对具有z值的数据进行子集，该数据匹配我想要的每个轮廓，并根据ID列查看它们与原始BW数据进行比较的纬度和经度。

Here is a picture of what I am talking about:

这是我正在谈论的图片：

I want to know which red points are within the 50% contour (blue) and which are within the 90% contour (red).

我想知道哪些红点在50％轮廓（蓝色）内，哪些在90％轮廓内（红色）。

Note: much of this code is from other questions. Big shout-out to all those who contributed!

注意：此代码的大部分来自其他问题。向所有贡献者致敬！

Thank you!

谢谢！

2 个解决方案

#1

You can use point.in.polygon from sp

你可以使用sp中的point.in.polygon

## Interactively check points
plot(bw)
identify(bw$lon, bw$lat, labels=paste("(", round(bw$lon,2), ",", round(bw$lat,2), ")"))

## Points within polygons
library(sp)
dens <- kde2d(x, y, n=200, lims=c(c(-73, -70), c(-13, -12)))  # don't clip the contour
ls <- contourLines(dens, level=levels)
inner <- point.in.polygon(bw$lon, bw$lat, ls[[2]]$x, ls[[2]]$y)
out <- point.in.polygon(bw$lon, bw$lat, ls[[1]]$x, ls[[1]]$y)

## Plot
bw$region <- factor(inner + out)
plot(lat ~ lon, col=region, data=bw, pch=15)
contour(dens, levels=levels, labels=prob, add=T)

#2

I think this is the best way I can think of. This uses a trick to convert the contour lines to SpatialLinesDataFrame objects using the ContourLines2SLDF() function from the maptools package. Then I use a trick outlined in Bivand, et al.'s Applied Spatial Data Analysis with R for converting the SpatialLinesDataFrame object to SpatialPolygons. These can then be used with the over() function to extract points within each contour polygon:

我认为这是我能想到的最佳方式。这使用一种技巧，使用maptools包中的ContourLines2SLDF（）函数将轮廓线转换为SpatialLinesDataFrame对象。然后我使用Bivand等人的应用空间数据分析与R中概述的技巧将SpatialLinesDataFrame对象转换为SpatialPolygons。然后可以将这些与over（）函数一起使用，以提取每个轮廓多边形内的点：

##  Simulate some lat/lon data:
x <- rnorm(363, 45, 10)
y <- rnorm(363, 45, 10)

##  Version 1 (without ggplot2):
library(MASS)
dens <- kde2d(x, y, n=200)

##  The contours to plot:
prob <- c(0.9, 0.5)
dx <- diff(dens$x[1:2])
dy <- diff(dens$y[1:2])
sz <- sort(dens$z)
c1 <- cumsum(sz) * dx * dy 
levels <- sapply(prob, function(x) { 
    approx(c1, sz, xout = 1 - x)$y
})
plot(x,y)
contour(dens, levels=levels, labels=prob, add=T)

##  Create spatial objects:
library(sp)
library(maptools)

pts <- SpatialPoints(cbind(x,y))

lines <- ContourLines2SLDF(contourLines(dens, levels=levels))

##  Convert SpatialLinesDataFrame to SpatialPolygons:
lns <- slot(lines, "lines")
polys <- SpatialPolygons( lapply(lns, function(x) {
  Polygons(list(Polygon(slot(slot(x, "Lines")[[1]], 
    "coords"))), ID=slot(x, "ID"))
    }))

##  Construct plot from your points, 
plot(pts)

##  Plot points within contours by using the over() function:
points(pts[!is.na( over(pts, polys[1]) )], col="red", pch=20)
points(pts[!is.na( over(pts, polys[2]) )], col="blue", pch=20)

contour(dens, levels=levels, labels=prob, add=T)