I'm using ArcGIS 10.4x and Python/numpy for this question. Can anyone suggest an efficient way to determine the array location of the nearest element whose value is less than the search element? As an example, if I provide the script a row/col number of a numpy array (the element value for example = 42.5), what is the most efficient way to find the row/col values for the nearest element whose value is less than the 42.5 value? I will have many source cells to iterate through and more often than not, an element with a lower value will be relatively close, so ideally I'd have a way where the search originates from the elements nearest to the source and radiates outwards until a solution is found.
2 Answers
You can do this with scipy (I hope it is included in ArcGis 10.4).
import numpy as np
from scipy.spatial.distance import cdist
a,b = np.where(array < 42.5) # get rows and columns where vlue is less than 42.5
x = zip(a,b) # create a list with (row,column)
d = np.argmin(cdist(np.array([[row,column]]), x)) # first distances are calculated between (row, col) of your input value, than nearest index value is selected
x[d] # this is result (row,col) of the nearest element under 42.5
-
Thanks very much for that. I'm be just had a peek and some machines are 10.3.1 (no scipy). Is there something comparable to scipy spatial.distance cdist() in numpy?– JaeCommented Aug 20, 2016 at 17:15
-
Also, is there an efficient way to limit the search distance besides creating a new array of a subset of the original? I.e. Only search within the adjacent 15 rows/Columns.– JaeCommented Aug 20, 2016 at 17:35
-
With numpy try:
np.linalg.norm(np.array(x)-np.array((row,col)), axis=1)– dmh126Commented Aug 20, 2016 at 17:37 -
And if you want to slice your 2D array, just do array[a:b, c:d]; ex. top left 3x4 array will be array[0:3,0:4]– dmh126Commented Aug 20, 2016 at 17:41
-
@Jae the more you use arcpy, the more you will want to also run regular python...I'd suggest having a separate install with all the goodies(but with no arcpy) for times like these. anaconda works well and you can install numpy/scipy/gdal super easily. Commented Aug 20, 2016 at 19:25
I converted point to raster, using field populated by 1 and tested:
arcpy.gp.Expand_sa("one", "D:/Scratch/two", "1", "1")
arcpy.gp.RasterCalculator_sa("""Con(IsNull("one"),"ORIGINAL" * "two")""", "D:/Scratch/three")
Originally it creates 8 neighbours around cell of interest.Repeat results in 16 cells etc.
This is rather simple thing to loop by assigning "ONE"="TWO" inside the script and using
arcpy.GetRasterProperties_management(in_raster="three", property_type="MINIMUM", band_index="")
to check if condition met and break. After that I suggest converting tiny raster "three" to points, delete ones that do not match criteria and do NEAR analysis.
