Jak zaimplementować dwuwymiarową funkcję K Ripleya?

9

Załączony obraz przedstawia szczelinę leśną z czerwoną sosną reprezentowaną przez koła i białą sosną reprezentowaną przez krzyże. Interesuje mnie ustalenie, czy istnieje pozytywny lub negatywny związek między dwoma gatunkami drzew sosny (tj. Czy rosną one na tych samych obszarach). Znam Kcrossa i Kmulti w pakiecie spatstat R. Ponieważ jednak mam 50 luk do przeanalizowania i jestem bardziej zaznajomiony z programowaniem w Pythonie niż R, chciałbym znaleźć iteracyjne podejście przy użyciu ArcGIS i Pythona. Jestem również otwarty na rozwiązania R.

Jak mogę zaimplementować dwuwymiarową funkcję K Ripleya?

wprowadź opis zdjęcia tutaj

Aaron
źródło
4
W drugim pytaniu możesz czerpać inspirację z tej odpowiedzi . W Pythonie tasowanie etykiet powinno być łatwe. W przypadku statystyk przestrzennych w Pythonie warto przyjrzeć się PySAL .
MannyG,

Odpowiedzi:

8

Po wielu poszukiwaniach w tylnych rogach dokumentacji ESRI doszedłem do wniosku, że nie ma rozsądnego sposobu uruchomienia dwuwymiarowej funkcji K Ripleya w Arcpy / ArcGIS. Znalazłem jednak rozwiązanie przy użyciu R:

# 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
źródło
3
Również FYI w bibliotece spatstat zastosowano dwuwariantową wersję Ripleya K. Niewłaściwe jest definiowanie obszaru badań za pomocą wypukłego kadłuba punktów, patrz funkcja ripras i cytowana literatura.
Andy W
2
Zauważ, że standaryzujesz oczekiwanie zerowe wokół zera, a tym samym wyprowadzasz statystyki Besag-L.
Jeffrey Evans
6

Jest wbudowane narzędzie skrypt o nazwie Multi-Odległość przestrzenna analiza skupień (Ripleys K Function) pod statystyka przestrzenna - Wzory Analizowanie ToolSet w ArcToolbox. Możesz przeczytać kod źródłowy narzędzia, jeśli przejdziesz do jego właściwości i zlokalizujesz skrypt używany na karcie Źródło.

blah238
źródło
Masz pomysł, jak uruchomić K jako funkcję dwuwymiarową w Arc, jeśli to w ogóle możliwe?
Aaron
1
Jestem pewien, że jest to możliwe, ale nie mogłem powiedzieć, jak to zrobić. Czy sprawdziłeś kod źródłowy wbudowanego narzędzia, aby zobaczyć, jakie zmiany należy wprowadzić?
blah238,
Kod źródłowy wygląda dość intensywnie. Zdecydowałem się zbadać rozwiązania R.
Aaron
3
Naprawdę nie zawracałbym sobie głowy próbą modyfikacji kodu ArcGIS Python. W najlepszym razie jest to kod spaghetti i nie wykonuje prawidłowego testu istotności. W przypadku problemów z procesem punktu dwuwymiarowego szczególnie ważne jest wykonanie testu istotności Monte Carlo, który jest dostępny w R z funkcją „obwiedni”.
Jeffrey Evans
1
Dzięki Jeffrey, nie wiem, o czym myślałem, polecając każdemu spojrzenie na kod źródłowy ESRI :)
blah238,