Como implementar a função K bivariada de Ripley?

9

A imagem em anexo mostra uma lacuna na floresta com pinheiro vermelho representado como círculos e pinheiro branco representado como cruzes. Estou interessado em determinar se existe uma associação positiva ou negativa entre as duas espécies de pinheiros (ou seja, se elas estão crescendo ou não nas mesmas áreas). Estou ciente do Kcross e do Kmulti no pacote Spstat R. No entanto, como tenho 50 lacunas para analisar e estou mais familiarizado com a programação em python do que R, gostaria de encontrar uma abordagem iterativa usando ArcGIS e python. Também estou aberto a soluções de R.

Como posso implementar a função K de Ripley bivariada?

insira a descrição da imagem aqui

Aaron
fonte
4
Para sua segunda consulta, você pode obter alguma inspiração dessa resposta . O embaralhamento de rótulos deve ser fácil no Python. Para estatísticas espaciais no Python, você pode querer dar uma olhada no PySAL .
precisa

Respostas:

8

Após muita pesquisa nos cantos traseiros da documentação da ESRI, concluí que não há uma maneira razoável de executar a função K de Ripley bivariada no Arcpy / ArcGIS. No entanto, eu encontrei uma solução 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
fonte
3
Além disso, para sua informação, a biblioteca spatstat possui uma implementação do K. bivariado de Ripley. É inadequado definir a área de estudo através do casco convexo dos pontos, ver a função ripras e a literatura citada.
Andy W
2
Observe que você está padronizando a expectativa nula em torno de zero e, portanto, derivando a estatística Besag-L.
Jeffrey Evans
6

Existe uma ferramenta de script embutida chamada Análise de Cluster Espacial Multi-Distância (Função Ripleys K) no conjunto de ferramentas Estatísticas Espaciais - Analisando Padrões no ArcToolbox. Você pode ler o código fonte da ferramenta se entrar em suas propriedades e localizar o script usado na guia Origem.

blah238
fonte
Alguma idéia de como executar K como uma função bivariada no Arc, se possível?
Aaron
11
Tenho certeza de que é possível, mas não sei dizer como fazê-lo. Você consultou o código-fonte da ferramenta interna para ver quais modificações precisam ser feitas?
blah238
O código fonte parece bastante intenso. Optei por explorar soluções R.
Aaron
3
Eu realmente não me incomodaria em tentar modificar o código do ArcGIS Python. É, na melhor das hipóteses, o código do espaguete e não realiza o teste de significância correto. Para problemas no processo de pontos bivariados, é especialmente importante executar um teste de significância de Monte Carlo, disponível em R com a função "envelop".
precisa saber é o seguinte
11
Graças Jeffrey, eu não sei o que eu estava pensando recomendando olhar ninguém no código fonte ESRI :)
blah238