Wie implementiere ich die bivariate Ripley K-Funktion?

9

Das beigefügte Bild zeigt eine Waldlücke mit roter Kiefer als Kreise und weißer Kiefer als Kreuze. Ich bin daran interessiert festzustellen, ob zwischen den beiden Kiefernarten eine positive oder negative Assoziation besteht (dh ob sie in denselben Gebieten wachsen oder nicht). Ich kenne Kcross und Kmulti im R-Spatstat-Paket. Da ich jedoch 50 Lücken zu analysieren habe und mit der Programmierung in Python besser vertraut bin als mit R, möchte ich einen iterativen Ansatz mit ArcGIS und Python finden. Ich bin auch offen für R-Lösungen.

Wie kann ich eine bivariate Ripley-K-Funktion implementieren?

Geben Sie hier die Bildbeschreibung ein

Aaron
quelle
4
Bei Ihrer zweiten Anfrage können Sie sich von dieser Antwort inspirieren lassen . Das Mischen von Etiketten sollte in Python einfach sein. Für räumliche Statistiken in Python sollten Sie sich PySAL ansehen .
MannyG

Antworten:

8

Nach langem Suchen in den hinteren Ecken der ESRI-Dokumentation bin ich zu dem Schluss gekommen, dass es keinen vernünftigen Weg gibt, eine bivariate Ripley-K-Funktion in Arcpy / ArcGIS auszuführen. Ich habe jedoch eine Lösung mit R gefunden:

# Calculates an estimate of the cross-type L-function for a multitype point pattern.
library(maptools)
library(spatstat)
library(sp)

# Subset certain areas within a points shapefile.  In this case, features are grouped by gap number
gap = 1

# Read the shapefile
sdata = readShapePoints("C:/temp/GapPoints.shp")  #Read the shapefile
data = sdata[sdata$SITE_ID == gap,]  # segregate only those points in the given cluster

# Get the convex hull of the study area measurements
gapdata = readShapePoints("C:/temp/GapAreaPoints_merged.shp")  #Read the shapefile that is used to estimate the study area boundary
data2 = gapdata[gapdata$FinalGap == gap,]  # segregate only those points in the given cluster
whole = coordinates(data2) # get just the coords, excluding other data
win = convexhull.xy(whole) # Convex hull is used to get the study area boundary
plot(win)

# Converting to PPP
points = coordinates(data) # get just the coords, excluding other data
ppp = as.ppp(points, win) # Convert the points into the spatstat format
ppp = setmarks(ppp, data$SPECIES) # Set the marks to species type YB or EH
summary(ppp) # General info about the created ppp object
plot(ppp) # Visually check the points and bounding area

# Plot the cross type L function
# Note that the red and green lines show the effects of different edge corrections
plot(Lcross(ppp,"EH","YB"))

# Use the Lcross function to test the spatial relationship between YB and EH
L <- envelope(ppp, Lcross, nsim = 999, i = "EH", j = "YB")
plot(L)
Aaron
quelle
3
Zu Ihrer Information, die Spatstat-Bibliothek hat eine Implementierung von bivariatem Ripley K. Es ist unangemessen, das Untersuchungsgebiet über die konvexe Hülle der Punkte zu definieren, siehe die Ripras- Funktion und die zitierte Literatur.
Andy W
2
Beachten Sie, dass Sie die Nullerwartung um Null standardisieren und somit die Besag-L-Statistik ableiten.
Jeffrey Evans
6

Unter dem Toolset " Räumliche Statistik - Analysieren von Mustern" in ArcToolbox befindet sich ein integriertes Skript-Tool namens " Räumliche Clusteranalyse mit mehreren Entfernungen" (Ripleys K-Funktion) . Sie können den Quellcode des Tools lesen, wenn Sie in dessen Eigenschaften das auf der Registerkarte Quelle verwendete Skript suchen.

blah238
quelle
Haben Sie eine Idee, wie Sie K als bivariate Funktion in Arc ausführen können, wenn dies überhaupt möglich ist?
Aaron
1
Ich bin mir sicher, dass es möglich ist, aber ich kann Ihnen nicht sagen, wie es geht. Haben Sie sich den Quellcode für das integrierte Tool angesehen, um festzustellen, welche Änderungen vorgenommen werden müssen?
blah238
Der Quellcode sieht ziemlich intensiv aus. Ich habe mich entschieden, R-Lösungen zu untersuchen.
Aaron
3
Ich würde wirklich nicht versuchen, den ArcGIS Python-Code zu ändern. Es ist bestenfalls Spaghetti-Code und führt nicht den richtigen Signifikanztest durch. Bei Problemen mit bivariaten Punktprozessen ist es besonders wichtig, einen Monte-Carlo-Signifikanztest durchzuführen, der in R mit der Funktion "Hüllkurve" verfügbar ist.
Jeffrey Evans
1
Danke Jeffrey, ich weiß nicht, was ich dachte, um jemandem zu empfehlen, sich den ESRI-Quellcode anzusehen :)
blah238