1

I’ve got a 2D NumPy Array derived from a Digital Elevation Model. I’m iterating through a series of specific row/col locations in the array and for each specified location, I’d like to define the nearest element whose value is lower than the current element I’m querying. Through the help of @dmh126 (Nearest numpy array element whose value is less than the current element), I’ve figured out the slice and the distance calculations, but I’m having trouble getting the coordinates of the “nearest” element from the sliced array back to the original. When the search element is near the center of the array, things work fine, but when the search element is near the axis edges, my slice entries become negative; which I believe can be a problem. I can’t seem to get the logic correct to transfer the coordinates of my “nearest” element from the sliced array back to the original for further processing. Any thoughts on what I’m doing wrong?

So far I’ve got:

# Window dimensions in metres
max_search_dist_m = 100
# Cell size of source DEM
CELLX = 1.0

# Array dimensions
arrayRowCount = array.shape[0]
arrayColCount = array.shape[1]

# Get the element value at the search location
theSinkZ = float(array[theSinkRow][theSinkCol])

# Define search window size
SinkWindow = int((float(max_search_dist_m) / float(CELLX)) / 2)

# Define corner elements for search window
col_min = theSinkCol - SinkWindow
col_max = theSinkCol + SinkWindow
row_min = theSinkRow - SinkWindow
row_max = theSinkRow + SinkWindow
col_min_orig = col_min
row_min_orig = row_min

# Define “windowed sink row/col” location. This is intended to be the row/col of the original search element relative to the sliced array.
if col_min < 0:
   col_min = 0
   WindowedSinkCol = SinkWindow + col_min_orig
else:
   WindowedSinkCol = SinkWindow
if col_max > arrayColCount: col_max = arrayColCount
if row_min < 0:
   row_min = 0
   WindowedSinkRow = SinkWindow + row_min_orig
else:
   WindowedSinkRow = SinkWindow

# Slice array to search window
# Get rows and columns where value is less than theSinkZ
c,e = np.where(array[int(col_min):int(col_max), int(row_min):int(row_max)] < theSinkZ) 

# Create a list of row/cols less than theSinkZ
y = zip(c,e)

# If solution exists, proceed.
if len(y) > 0:
   # first distances are calculated between (row, col) of your input value, than nearest index value is selected.
   d = np.argmin(cdist(np.array([[WindowedSinkRow,WindowedSinkCol]]), y))

   # Extract nearest row/col location of Sliced array
   WindowedOutCol = y[d][1]
   WindowedOutRow = y[d][0]

   # Identify outlet location relative to original full array
   theOutCol = theSinkCol + col_min
   theOutRow = theSinkRow + row_min
3

It is easy if you understand the slicing in numpy

import numpy as np
x = range(16)
x = np.reshape(x,(4,4))
print x
[[ 0  1  2  3]
[ 4  5  6  7]
[ 8  9 10 11]
[12 13 14 15]]
b = x[2:4,1:4]
print b
[[ 9 10 11]
[13 14 15]]

Now, in an easiest way:

print  b[0][0]
9
# is equivalent to
print x[0+2][0+1]
9

Therefore simplifying your script

# conversion sliced to array indices
def col(i):
   return i+col_min
def row(i):
   return i+row_min
from osgeo import gdal
g = gdal.Open("a.tif")
array = g.ReadAsArray()
theSinkZ = array[30][50]
arrayRowCount, arrayColCount= array.shape
theSinkRow, theSinkCol = 30,50
SinkWindow = 5
# Define corner elements for search window
col_min = theSinkCol - SinkWindow
col_max = theSinkCol + SinkWindow
row_min = theSinkRow - SinkWindow
row_max = theSinkRow + SinkWindow

enter image description here

y = zip(*np.where(array[int(col_min):int(col_max), int(row_min):int(row_max)]  > 220))

The resulting points are in green

enter image description here

from scipy.spatial.distance import cdist
WindowedSinkCol = SinkWindow
WindowedSinkRow = SinkWindow
d = np.argmin(cdist(np.array([[WindowedSinkRow,WindowedSinkCol]]), y))
a, b = y[d]
# sliced indices
print a, b
9 13
print array[int(col_min):int(col_max), int(row_min):int(row_max)][a][b]
220.07
# array indices
print col(a), row(b)
49 33
print array[col(a)][row(b)]
220.07

Result: nearest point (in orange)

enter image description here

But there are other solutions as using masks (with Numpy or GDAL) for example.

New

The slicing use only the values in the original array (no need of if...)

x = [[ 0  1  2  3]
[ 4  5  6  7]
[ 8  9 10 11]
[12 13 14 15]]
x.shape
(4,4)
b = x[2:4,1:4]
print b
[[ 9 10 11]
[13 14 15]]
b = x[2:8,1:4]
print b
[[ 9 10 11]
[13 14 15]]
b = x[-2:8,1:8]
print b
[[ 9 10 11]
[13 14 15]]

A confirmation is given if I save the sliced array (array[row_min:rox_max, col_min:col_max], Numpy use (row, col) ordering.) as a new raster

enter image description here

And if I apply your script

enter image description here

  • Thanks @gene for this response. When I'm processing a row/column location that is near the origin, say at [2,3], the sliced array (assuming SinkWindow = 5) would have a row_min of 2-5 = -3 and col_min of 3-5 = -2. Can I ignore this? I had assumed that I would require an if statement to ensure I was within bounds of the array. And if I need to set these values to 0, WindowedSinkRow and WindowedSinkCol would not equal SinkWindow. Right? – Jae Sep 4 '16 at 19:25
  • look above in new – gene Sep 5 '16 at 15:46
  • Thanks @gene. It's great that python can handle "out of bounds" slice ranges, but what I'm getting at is if your last example, the "focus cell" would no longer be equal to "SinkWindow" since the sliced array is not equally surrounding the "focus cell". So it seems to me that if row or col min < 0, wouldn't I just use WindowedSinkRow = SinkRow instead of WindowedSinkRow = SinkWindow? P.S. I wish I was in front of a computer so I could test instead of asking silly questions. Just trying to make sure I understand the logic properly. – Jae Sep 5 '16 at 18:43

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.