Comment implémenter la fonction K bivariée de Ripley?

9

L'image ci-jointe montre un espace forestier avec du pin rouge représenté comme des cercles et du pin blanc représenté comme des croix. Je souhaite déterminer s'il existe une association positive ou négative entre les deux espèces de pins (c'est-à-dire si elles poussent ou non dans les mêmes zones). Je connais Kcross et Kmulti dans le package R spatstat. Cependant, comme j'ai 50 lacunes à analyser et que je connais mieux la programmation en python que R, je voudrais trouver une approche itérative utilisant ArcGIS et python. Je suis également ouvert aux solutions R.

Comment puis-je implémenter une fonction K bivariée de Ripley?

entrez la description de l'image ici

Aaron
la source
4
Pour votre deuxième demande, vous pourriez vous inspirer de cette réponse . Le brassage des étiquettes devrait être facile en Python. Pour les statistiques spatiales en Python, vous voudrez peut-être regarder PySAL .
MannyG

Réponses:

8

Après de nombreuses recherches dans les coins arrière de la documentation ESRI, j'ai conclu qu'il n'y avait aucun moyen raisonnable d'exécuter une fonction K bivariée de Ripley dans Arcpy / ArcGIS. Cependant, j'ai trouvé une solution en utilisant 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
la source
3
Pour info, la bibliothèque spatstat a une implémentation du K. bivarié de Ripley. Il est inapproprié de définir la zone d'étude via la coque convexe des points, voir la fonction ripras et la littérature citée.
Andy W
2
Notez que vous standardisez l'espérance nulle autour de zéro et dérivez ainsi la statistique Besag-L.
Jeffrey Evans
6

Il existe un outil de script intégré appelé Analyse de cluster spatial à distances multiples (fonction Ripleys K) sous le jeu d'outils Statistiques spatiales - Analyse des modèles dans ArcToolbox. Vous pouvez lire le code source de l'outil si vous allez dans ses propriétés et localisez le script utilisé dans l'onglet Source.

blah238
la source
Une idée sur la façon d'exécuter K en tant que fonction bivariée dans Arc, si possible?
Aaron
1
Je suis sûr que c'est possible, je ne pourrais pas vous dire comment le faire cependant. Avez-vous regardé le code source de l'outil intégré pour voir quelles modifications doivent être apportées?
blah238
Le code source semble assez intense. J'ai choisi d'explorer les solutions R.
Aaron
3
Je ne prendrais vraiment pas la peine d'essayer de modifier le code ArcGIS Python. Il s'agit au mieux de code spaghetti et n'effectue pas le test de signification correct. Pour les problèmes de processus ponctuels bivariés, il est particulièrement important d'effectuer un test de signification Monte Carlo, qui est disponible en R avec la fonction «enveloppe».
Jeffrey Evans
1
Merci Jeffrey, je ne sais pas ce que je pensais recommander à quiconque de regarder le code source ESRI :)
blah238