Transforming Indices from NumPy Array to Sliced Array, then Back

Question

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

Answer 1

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

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

The resulting points are in green

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)

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

And if I apply your script