आर - विशिष्ट कंटूर

Question

के भीतर अंक कैसे प्राप्त करें मैं लेट और लॉन डेटा पर kde2d (MASS) के साथ घनत्व प्लॉट बना रहा हूं। मैं जानना चाहता हूं कि मूल डेटा से कौन से बिंदु एक विशिष्ट समोच्च के भीतर हैं।

मैं दो दृष्टिकोणों का उपयोग करके 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) 

और यहां 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") 

मैं प्लॉट और लोन पॉइंट प्लॉट और 90% और 50% समोच्चों के साथ भूखंड बनाते हैं। मैं बस जानना चाहता हूं कि 90% और 50% समोच्चों के भीतर सटीक बिंदुओं को कैसे निकालना है।

मैंने z मानों (kde2d से घनत्व भूखंडों की ऊंचाई) को खोजने की कोशिश की है जो लेट और लोन मानों की प्रत्येक पंक्ति से जुड़े हुए हैं लेकिन कोई भाग्य नहीं था। मैं यह भी सोच रहा था कि मैं प्रत्येक पंक्ति को लेबल करने के लिए डेटा में एक आईडी कॉलम जोड़ सकता हूं और फिर melt() का उपयोग करने के बाद इसे किसी भी तरह स्थानांतरित कर सकता हूं। फिर मैं उस डेटा को बस सब्सक्राइब कर सकता हूं जिसमें ज़ेड के मान हैं जो प्रत्येक समोच्च से मेल खाते हैं और देखें कि कौन सी लेट और लोन की तुलना आईडी कॉलम के आधार पर मूल बीडब्ल्यू डेटा से की जाती है।

यहाँ है कि मैं क्या बात कर रहा हूँ की एक तस्वीर है:

मुझे पता है जो लाल अंक 50% समोच्च (नीला) के भीतर कर रहे हैं और जो 90% समोच्च भीतर (लाल कर रहे हैं)।

नोट: यह कोड अन्य प्रश्नों से है। योगदान देने वालों के लिए बड़ा चिल्लाओ!

धन्यवाद!

Answer 1

आप से 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

मुझे लगता है कि यह सबसे अच्छा तरीका है जिसके बारे में मैं सोच सकता हूं। यह maptools पैकेज से ContourLines2SLDF() फ़ंक्शन का उपयोग करके समोच्च रेखाओं को SpatialLinesDataFrame ऑब्जेक्ट्स में कनवर्ट करने के लिए एक चाल का उपयोग करता है। फिर मैं SpatialLinesDataFrame ऑब्जेक्ट को SpatialPolygons में कनवर्ट करने के लिए आर के साथ Bivand, et al। एप्लाइड स्पेटियल डेटा विश्लेषण में उल्लिखित एक चाल का उपयोग करता हूं। ये तो 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)