R - Wie finde ich Punkte innerhalb einer bestimmten Kontur?

Question

Ich erstelle Dichtendiagramme mit kde2d (MASS) für lat und lon Daten. Ich würde gerne wissen, welche Punkte aus den Originaldaten innerhalb einer bestimmten Kontur liegen.

Ich erstellen 90% und 50% Konturen mit zwei Ansätze. Ich möchte wissen, welche Punkte innerhalb der 90% -Kontur liegen und welche Punkte innerhalb der 50% -Kontur liegen. Die Punkte in der 90% -Kontur enthalten alle innerhalb der 50% -Kontur. Der letzte Schritt besteht darin, die Punkte innerhalb der 90% -Kontur zu finden, die nicht innerhalb der 50% -Kontur liegen (ich brauche bei diesem Schritt nicht unbedingt Hilfe).

# 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) 

Und hier ist Version 2 - mit ggplot2. Ich würde diese Version am liebsten verwenden, um die Punkte innerhalb der Konturen 90% und 50% zu finden.

# 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") 

Ich erstelle die Plots mit allen lat und lon Punkten gezeichnet und 90% und 50% Konturen. Ich möchte einfach wissen, wie man die genauen Punkte extrahiert, die innerhalb der 90% und 50% Konturen liegen.

Ich habe versucht, die z-Werte (die Höhe der Dichteplots von kde2d) zu finden, die mit jeder Reihe von Lat- und lon-Werten verbunden sind, aber kein Glück hatten. Ich dachte auch, dass ich eine ID-Spalte zu den Daten hinzufügen könnte, um jede Zeile zu beschriften und diese dann nach der Verwendung von melt() irgendwie zu übertragen. Dann könnte ich einfach die Daten, die Werte von z haben, die jeder gewünschten Kontur entsprechen, unterteilen und sehen, welches Breiten- und Längenverhältnis sie mit den ursprünglichen BW-Daten basierend auf der ID-Spalte vergleichen. Hier

ist ein Bild von dem, was ich rede:

Ich mag wissen, welche roten Punkte sind innerhalb der 50% Kontur (blau) und die innerhalb der 90% Kontur (rot).

Hinweis: ein Großteil dieses Codes stammt aus anderen Fragen. Großes Dankeschön an alle, die dazu beigetragen haben!

Vielen Dank!

Answer 1

Sie können point.in.polygon von sp

## 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) 

Answer 2

Ich denke, das ist der beste Weg, den ich mir vorstellen kann. Dies verwendet einen Trick, um die Konturlinien in SpatialLinesDataFrame Objekte mit der ContourLines2SLDF() Funktion aus dem maptools Paket zu konvertieren. Dann verwende ich einen Trick, der in Bivand, et al., Angewandte räumliche Datenanalyse mit R zum Umwandeln des SpatialLinesDataFrame Objekts zu SpatialPolygons umrissen ist. Diese können dann mit der over() Funktion verwendet werden, um Punkte innerhalb jeder Kontur Polygon zu extrahieren:

## 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)