# How to implement bivariate Ripley's K function?

The attached image shows a forest gap with red pine represented as circles and white pine represented as crosses. I am interested in determining if there is a positive or negative association between the two species of pine trees (i.e. whether or not they are growing in the same areas). I am aware of Kcross and Kmulti in the R spatstat package. However, since I have 50 gaps to analyze and am more familiar with programming in python than R, I would like to find an iterative approach using ArcGIS and python. I am also open to R solutions.

How can I implement a bivariate Ripley's K function? • For your second inquiry, you might glean some inspiration from this answer. The shuffling of labels should be easy in Python. For spatial stats in Python you might want to look at PySAL. – MannyG Dec 29 '12 at 23:42

After much searching in the back corners of ESRI documentation, I've concluded that there is no reasonable way run a bivariate Ripley's K function in Arcpy/ArcGIS. However, I have found a solution using 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

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)
``````
• Also FYI the spatstat library has an implementation of bivariate Ripley's K. It is inappropriate to define the study area via the convex hull of the points, see the ripras function and the cited literature. – Andy W Dec 31 '12 at 1:57
• Note that you are standardizing the null expectation around zero and thus deriving the Besag-L statistic. – Jeffrey Evans Jan 31 '13 at 18:12

There is a built-in script tool called Multi-Distance Spatial Cluster Analysis (Ripleys K Function) under the Spatial Statistics - Analyzing Patterns toolset in ArcToolbox. You can read the tool's source code if you go into its properties and locate the script used in the Source tab.

• Any idea on how to run K as a bivariate function in Arc, if at all possible? – Aaron Dec 29 '12 at 20:37
• I am sure it is possible, I couldn't tell you how to do it though. Have you looked at the source code for the built-in tool to see what modifications need to be made? – blah238 Dec 29 '12 at 22:04
• The source code looks pretty intense. I've opted to explore R solutions. – Aaron Dec 30 '12 at 19:36
• I would really not bother trying to modify the ArcGIS Python code. It is spaghetti code at best and does not perform the correct significance test. For bivariate point process problems, it is especially important to perform a Monte Carlo significance test, which is available in R with the "envelop" function. – Jeffrey Evans Jan 31 '13 at 18:10
• Thanks Jeffrey, I don't know what I was thinking recommending anyone look at ESRI source code :) – blah238 Jan 31 '13 at 18:11