R - Come trovare punti all'interno di Contour specifico

Question

Sto creando grafici di densità con kde2d (MASS) su dati lat e lon. Vorrei sapere quali punti dei dati originali si trovano all'interno di un contorno specifico.

Creo i contorni del 90% e 50% utilizzando due approcci. Voglio sapere quali sono i punti all'interno del contorno del 90% e quali sono all'interno del contorno del 50%. I punti nel contorno del 90% conterranno tutti quelli all'interno del contorno del 50%. Il passo finale è trovare i punti all'interno del contorno del 90% che non si trovano all'interno del contorno del 50% (non ho necessariamente bisogno di aiuto con questo passaggio).

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

E qui è la versione 2 - utilizzando ggplot2. Preferirei utilizzare questa versione per trovare i punti all'interno dei contorni del 90% e 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") 

creo le trame con tutti i punti lat e lon tracciati e il 90% e il 50% contorni. Voglio semplicemente sapere come estrarre i punti esatti all'interno dei contorni del 90% e 50%.

Ho provato a trovare i valori z (l'elevazione dei grafici di densità da kde2d) associati a ciascuna riga di valori lat e lon ma non ha avuto fortuna. Stavo anche pensando di aggiungere una colonna ID ai dati per etichettare ciascuna riga e poi trasferirla in qualche modo dopo aver usato melt(). Quindi potrei semplicemente impostare sottoinsieme i dati che hanno valori di z che corrispondono a ciascun contorno che voglio e vedere quale lat e lon sono confrontati con i dati BW originali basati sulla colonna ID.

Ecco una foto di quello di cui sto parlando:

Voglio sapere che i punti rossi sono all'interno del 50% di contorno (blu) e quali sono all'interno del contorno del 90% (rosso).

Nota: gran parte di questo codice proviene da altre domande. Grande grido a tutti quelli che hanno contribuito!

Grazie!

Answer 1

È possibile utilizzare point.in.polygon da 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

Penso che questo sia il modo migliore che riesco a pensare. Questo utilizza un trucco per convertire le curve di livello in oggetti SpatialLinesDataFrame usando la funzione ContourLines2SLDF() dal pacchetto maptools. Quindi utilizzo un trucco delineato nell'analisi dei dati spaziali applicati di Bivand, et al. con R per la conversione dell'oggetto SpatialLinesDataFrame in SpatialPolygons. Questi possono poi essere utilizzati con la funzione over() per estrarre i punti all'interno di ogni poligono di contorno:

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