0

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.

1

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? – Jae Aug 20 '16 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. – Jae Aug 20 '16 at 17:35
  • With numpy try: np.linalg.norm(np.array(x)-np.array((row,col)), axis=1) – dmh126 Aug 20 '16 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] – dmh126 Aug 20 '16 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. – user1269942 Aug 20 '16 at 19:25
1

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. enter image description here

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.