¿Cómo implementar la función bivariada K de Ripley?

9

La imagen adjunta muestra una brecha forestal con pino rojo representado como círculos y pino blanco representado como cruces. Estoy interesado en determinar si existe una asociación positiva o negativa entre las dos especies de pinos (es decir, si crecen o no en las mismas áreas). Estoy al tanto de Kcross y Kmulti en el paquete R spatstat. Sin embargo, dado que tengo 50 lagunas para analizar y estoy más familiarizado con la programación en python que con R, me gustaría encontrar un enfoque iterativo usando ArcGIS y python. También estoy abierto a soluciones R.

¿Cómo puedo implementar una función K de Ripley bivariada?

ingrese la descripción de la imagen aquí

Aaron
fuente
44
Para su segunda consulta, puede obtener algo de inspiración de esta respuesta . La mezcla de etiquetas debería ser fácil en Python. Para las estadísticas espaciales en Python, es posible que desee mirar PySAL .
MannyG

Respuestas:

8

Después de mucho buscar en las esquinas traseras de la documentación de ESRI, he concluido que no hay una forma razonable de ejecutar una función K de Ripley bivariada en Arcpy / ArcGIS. Sin embargo, he encontrado una solución usando 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
fuente
3
También para su información, la biblioteca de Sstastat tiene una implementación de Ripley's K. bivariada. Es inapropiado definir el área de estudio a través del casco convexo de los puntos, ver la función de Ripras y la literatura citada.
Andy W
2
Tenga en cuenta que está estandarizando la expectativa nula en torno a cero y, por lo tanto, deriva la estadística Besag-L.
Jeffrey Evans el
6

Existe una herramienta de secuencia de comandos incorporada llamada Análisis de conglomerados espaciales de distancias múltiples (función Ripleys K) en el conjunto de herramientas Estadísticas espaciales - Análisis de patrones en ArcToolbox. Puede leer el código fuente de la herramienta si accede a sus propiedades y localiza el script utilizado en la pestaña Fuente.

blah238
fuente
¿Alguna idea sobre cómo ejecutar K como una función bivariada en Arc, si es posible?
Aaron
1
Estoy seguro de que es posible, aunque no podría decirte cómo hacerlo. ¿Has mirado el código fuente de la herramienta incorporada para ver qué modificaciones deben hacerse?
blah238
El código fuente se ve bastante intenso. He optado por explorar las soluciones R.
Aaron
3
Realmente no me molestaría en tratar de modificar el código Python de ArcGIS. Es un código de espagueti en el mejor de los casos y no realiza la prueba de significación correcta. Para problemas de proceso de puntos bivariados, es especialmente importante realizar una prueba de significancia de Monte Carlo, que está disponible en R con la función "envoltura".
Jeffrey Evans el
1
Gracias Jeffrey, no sé lo que estaba pensando en recomendar que alguien mire el código fuente de ESRI :)
blah238