Hi Dan Patterson,
I am at a stopping point for now as I did not have much time to work on this as it is not a "have to" since I currently have a tool that calls an excellent and fully debugged R-script (R code written by the original developer of the model) to do the job. But, it would be nice to be able do it in Arc to avoid an Arc plus an R install (although the R version with R script version install is not a big deal since multiple versions of R can be installed). My first ArcGIS version that uses the PRO kde kernel does a pretty good job with the arcpy kernel (KDE) and creates a nice polygon (with the use of some additional tools), but it would be nice to be able to use the gaussian kernel. This test routine does create a grid but it is not correct so it is only a start. What I have learned is that the documentation for arcpy is excellent but the documentation for SCIPY is not so great (at least from my learning perspective). So, this will have to sit for a while as I have some new work to get done thus I am not able to close this discussion with a solution.
DRAFT CODE SO FAR FOLLOWS:
def kerneldensityhrd(homerangedir, HomeRangeGDB, Points, SelBand, SeaRadUse, Cellsize, DebugThis, HrefMult):
if SelBand == "Reference" or SelBand == "Reference2":
# Determine bandwidth for Href (Reference)
arcpy.AddXY_management(Points)
xlist = []
ylist = []
href_curs = arcpy.da.SearchCursor(Points, ["POINT_X","POINT_Y"])
for row in href_curs:
pointx = row[0]
pointy = row[1]
xlist.append(pointx)
ylist.append(pointy)
Href = Hreference(xlist,ylist,DebugThis)
if DebugThis == "ON":
arcpy.AddMessage("Href = " + str(Href))
searchRadius = (Href + 2000) * HrefMult
else: #ArcGIS and ArcGIS2
searchRadius = SeaRadUse
# Execute KernelDensity
#Gaussian Kernel Density from point layer
KDEtest = "Yes"
if KDEtest == "Yes":
cellSize = int(Cellsize)
#Create paired coordinate list from point layer
my_list1 = []
my_list2 = []
select_curs = arcpy.da.SearchCursor(Points, ["POINT_X","POINT_Y"])
for row in select_curs:
xcoord = float(row[0])
ycoord = float(row[1])
#xycoord = '[ ' + str(xcoord) + ', ' + str(ycoord) + ']'
my_list1.append(xcoord)
my_list2.append(ycoord)
#arcpy.AddMessage(my_list1)
#arcpy.AddMessage(my_list2)
ymin = min(my_list2)
ymax = max(my_list2)
xmin = min(my_list1)
xmax = max(my_list1)
#temp = list(zip(my_list1,my_list2))
#my_array = np.array(temp)
m1 = my_list1
m2 = my_list2
X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([m1, m2])
kernel = stats.gaussian_kde(values, bw_method=searchRadius)
Z = np.reshape(kernel(positions).T, X.shape)
outKernelDensity = arcpy.NumPyArrayToRaster(Z, arcpy.Point(xmin,ymin), cellSize)
#GRAPH
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
extent=[xmin, xmax, ymin, ymax])
ax.plot(m1, m2, 'k.', markersize=2)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
#plt.show()
#exit(0)
else:
