topic Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange in ArcMap Questions

ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Tue, 21 Aug 2018 14:29:52 GMT

Hello,

A kernel density routine developed for the "R" language has been popular for use by government agencies to estimate core home range for Bighorn Sheep using telemetry points as input data. Mostly out of curiosity, I developed a routine using the most recent ArcGIS kernel density routine in Spatial Analyst. It produces very similar results using the HREF bandwidth. I asked the developer of the R home range routine what he thought and this was his answer:

This is his reasoning for using only the R routine and kernel density function and not using the ArcGIS python routine that I developed:

"The main issue is that the Home Range Arc tool uses a quartic kernel density function which is only an approximation to the Gaussian density function employed by the R home range estimator. ArcGIS uses it because it's fast and good enough for visualizations for every day users, but it's not a density function that you actually see in the home range literature. It's also not the density function used when we were estimating the foray frequencies and distances, which is in itself a reason not to go with it."

So, my question to you is "do you agree that ArcGIS is not capable of generating a home range polygon from telemetry points?". My reasoning is that it is because coming up with a representative polygon is an art and not an exact science. Sometimes the home range polygon is even hand delineated using expert knowledge. Also, the ArcGIS kernel density function was improved with version 10.3.

When using the HREF bandwidth what I have observed is that the R-script performance is faster on small datasets of telemetry points but the ArcGIS python routine is faster on larger datasets. The resulting polygon is very close.

Your thoughts regarding this will be appreciated. I have the routine working with both desktop and PRO.

Thanks,

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Tue, 21 Aug 2018 17:54:25 GMT

everyone has an opinion.

you can roll out your own if you are concerned since you have access to numpy and scipy in arcmap/ArcGIS pro

scipy.stats.gaussian_kde — SciPy v1.1.0 Reference Guide

some of the same literature is referenced and you could do a comparison if you wanted

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Wed, 22 Aug 2018 03:22:43 GMT

Dan,

Thanks for the good information. I will look into these new

capabilities.

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Sat, 11 Dec 2021 12:41:26 GMT

you can try their data for starters... just to give you an idea how to set it up with some 'fake' parameters

from scipy import stats

def measure(n):
    "Measurement model, return two coupled measurements."
    m1 = np.random.normal(size=n)
    m2 = np.random.normal(scale=0.5, size=n)
    return m1+m2, m1-m2


n = 5

m1, m = measure(n)

m1, m2 = measure(n)

xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()



X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]

positions = np.vstack([X.ravel(), Y.ravel()])

values = np.vstack([m1, m2])

kernel1 = stats.gaussian_kde(values, 'scott')

kernel2 = stats.gaussian_kde(values, 'silverman')

Z1 = np.reshape(kernel1(positions).T, X.shape)

Z2 = np.reshape(kernel2(positions).T, X.shape)

np.allclose(Z1, Z2)

True  # result‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Thu, 23 Aug 2018 19:18:55 GMT

Hi All,

I worded this sentence incorrectly; "do you agree that ArcGIS is not capable of generating a home range polygon from telemetry points?" I should have said; "do your agree that the ArcGIS 10.5.1 kernel ( quartic kernel ) should not be used for generating a bighorn sheep home range polygon from telemetry points since a true Gaussian Kernel is preferred".

And, Dan Patterson, it looks like scipy may be the solution to a true Gaussian kernel from ArcGIS Pro assuming that I can return the result as a raster registered to the real world. My existing routine uses only arcpy and spatial analyst plus calculated Href for bandwidth (has other bandwidth options) and does not use either scipy or numpy. It currently works with desktop versions 10.3.1 up to 10.6. It works with earlier versions of ArcGIS PRO but has a problem with a missing OID in PRO 2.2.1. I use several spatial analyst tools and some seem to work inconsistently with different versions of ArcGIS so I have found solutions that work for all of them (except now, PRO 2.2.1).

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Thu, 23 Aug 2018 20:25:57 GMT

creating the raster with coordinates is no problem

NumPyArrayToRaster—ArcPy Functions | ArcGIS Desktop

check my blog through my profile... I have tons of numpy (scipy etc) posts

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Fri, 31 Aug 2018 17:43:14 GMT

Hi Dan P.,

Thanks for your examples.

I have an ArcGIS raster to numpy array and back to ArcGIS raster working via the following code:

from scipy import stats
import numpy as np
import arcpy

#Input and output parameters
inRas = arcpy.GetParameterAsText(0)
outFolder = arcpy.GetParameterAsText(1)

# Set environment
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = inRas
arcpy.env.cellSize = inRas

# Lower left coordinate of block (in map units)
inRasDesc = arcpy.Describe(inRas)
result_left = arcpy.GetRasterProperties_management(inRas, "LEFT")
result_left = float(str(result_left))
result_bot = arcpy.GetRasterProperties_management(inRas, "BOTTOM")
result_bot = float(str(result_bot))
cellSize = inRasDesc.MeanCellWidth

my_array = arcpy.RasterToNumPyArray(inRas)
outRas = arcpy.NumPyArrayToRaster(my_array, arcpy.Point(result_left,result_bot), cellSize)
outRas.save(outFolder + "\\TestOut")

The output result exactly aligns with the input raster.

I also have a point feature class to numpy array working as prework to use the gaussian KDE density based on telemetry points.

my_array = arcpy.da.FeatureClassToNumPyArray (Location + "\\\\" + "Z" + str(animal) + "pts", ["OID@", "SHAPE@XY"])

However, the documentation says you can use a scalar. So, I'm creating a reference bandwidth (rather than using the default) in a separate routine and trying to input that as a scalar.

#Href formula is
# The square root ((varx + vary) / 2) times n to the -1/6 power

. . . etc,

return Href

However I don't understand either the syntax or the scalar data type. searchRadius is passed as a floating point number.

I have tried: kernel = stats.gaussian_kde(my_array, bw_method=searchRadius) and other variations but still get errors.

The searchRadius is a multiplier of my computed Href

Your thoughts?

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Fri, 31 Aug 2018 20:23:46 GMT

quick look for now

inRas must be converted to a float because it is text right now

/blogs/dan_patterson/2016/08/14/script-formatting Dan can you reformat the code so it has line numbers etc.

Parameters will be text using GetParameterAsText and if they are supposed to be int or float, they need to be converted appropriately

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Fri, 31 Aug 2018 21:59:44 GMT

Hi Dan,

I am sorry for the confusion. That first script using inRas has nothing to do with my problem. It works well. The arguments are only used to locate the raster and the output GDB to save to. I was just giving you some feedback that I was successful in being able to go from a geodatabase raster to a numpy array and then back to an identical geodatabase raster per my question in another post. I.E. that first script works well for me solving a problem that I asked you about in another post below this one (creating a raster with coordinates).

The 2nd part of the post above is the problem. I'm trying to pass a known bandwidth to the scipy gaussian kernel density routine rather than using one of the other options. I do this with the ArcGIS kernel density routine by specifying a calculated searchRadius. However, the gaussian routine comes up with an error on the bandwidth option when I use:

kernel = stats.gaussian_kde(my_array, bw_method=searchRadius)

Search radius is a calculated floating point number using the reference bandwidth equation. It could be that I am fooled by the bandwidth parameter error message and that the array I am passing is messing it up. The array is generated from a telemetry point feature class or shape file (a selection of points for a particular animal) and the following statement seems to work to successfully create the numpy array:

my_array = arcpy.da.FeatureClassToNumPyArray (Location + "\\\\" + "Z" + str(animal) + "pts", ["OID@", "SHAPE@XY"])

But, then when passing the array and search radius to the gaussian kde that statement fails saying there is a problem with the bandwidth parameter (bw_estimator).

There is no hurry on my end and you don't need to worry about this - just thought I would throw out the problem in case anyone knew a quick answer.

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Sat, 01 Sep 2018 00:21:45 GMT

Ok... that array is not of the form that is required by the kernel density estimate.

The current array is called a 'structured array' and since it has the id and the X,Y coordinates, there is a little trick to extract the x,y coordinates so that you get an 'ndarray' which is the floating point values of the array without the field names.

I have several little functions that make it easier, and I can pop them into a documentation example if you zip a shapefile with sample of the data so I can document the steps to make it easier.

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Sat, 01 Sep 2018 02:09:47 GMT

Hi Dan P.,

The data is sensitive so I cannot share it. It is just a regular point feature class with some attributes. But, you have me on the right track now. I will double check and use a non scaler (default without bandwidth) and it should flag the array as the error or even still the bandwidth even though only one parameter. Then I will look into array types or perhaps make an numpy array from my list of x,y points that I use to get the reference value. So, it must think the 3rd value in the existing array is band width parameter. I just need to put a little more effort into finding a solution.

Thanks,

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Sat, 01 Sep 2018 03:18:57 GMT

Ok... I will send some fake point data so you can see how it is extracted

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Mon, 03 Sep 2018 18:54:31 GMT

Hey Dan P.,

Thanks, but no need to send the sample data. The clue you gave me was sufficient for me to get past that step and move on to the next one.

Here is the rough routine of getting Gaussian kernel density from a point layer (in this case a shape file). Note that the searchRadius variable mentioned here was passed from another routine that calculated it using the reference bandwidth.

      #Gaussian Kernel Density from point layer
      import arcpy
      from scipy import stats
      import numpy as np

      KDEtest = "Yes"
      if KDEtest == "Yes":
         #Create paired coordinate list from point layer
         my_list1 = []
         my_list2 = []
         select_curs = arcpy.da.SearchCursor(Location + "\\\\" + "Z" + str(animal) + "pts", ["POINT_X","POINT_Y"], query)
         for row in select_curs: #Cursor is an extract from point layer for an individual animal
            xcoord = float(row[0])
            ycoord = float(row[1])
            my_list1.append(xcoord)
            my_list2.append(ycoord)
         temp = list(zip(my_list1,my_list2)) #The zip pairs the x,y coordinates into one list

#Check list result
arcpy.AddMessage(temp)

         #Create Numpy array from list and check result
         my_array = np.array(temp)
         arcpy.AddMessage(my_array)

#Gaussian Kernel Density routine
kernel = stats.gaussian_kde(my_array, bw_method=searchRadius)

#Next challenge is to figure out how to turn the kde kernel into a grid?

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Tue, 04 Sep 2018 22:45:49 GMT

Dan

sounds good

Post your solution back when you get done and mark this closed or change to a Discussion.

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Tue, 02 Oct 2018 14:56:38 GMT

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:

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Tue, 02 Oct 2018 15:02:44 GMT

Dan... I will have a look some time soon. And yes, SciPy's documentation could use some work.

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanGreene — Thu, 14 Feb 2019 16:07:56 GMT

Hi Dan P.,

You could close this topic if you would like to since I don't have time to dig into it anymore. The ideal would be if ESRI supported the standard kernel in addition to the existing one. Interest seems to have died down and calling an R-script does the job for now.

Thanks,

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanielGreene — Wed, 03 Jul 2019 01:33:36 GMT

Thanks Dan P.,

I don't know if I will have time down the road to get back to this. My thoughts are that my python routine using the existing ArcGIS kernel to identify a home range polygon works well and is close to the one generated via the desired kernel. In reality it is only the scientists who object. But in reality, there is no absolute correct home range polygon and it could just as well be delineated by a wildlife biologist. Who is really to say what is the best polygon but I am not a kernel density expert.

Dan

Re: ArcGIS 10.5.1 Kernel Density and Bighorn Sheep Homerange

DanPatterson_Retired — Wed, 03 Jul 2019 01:43:15 GMT

So true...