How to identify the highest and lowest pixels within a defined neighborhood area of a central pixel
So, let me get this straight, for every 3x3 moving window, you want the location of the min within that moving window?
Lets look at the top 6 rows and columns. Here they are in array form. Array 'a' is the original dem, with the top row being a nodata row)
a[:6, :6] # original dem with nodata row
array([[-32768, -32768, -32768, -32768, -32768, -32768],
[ 619, 606, 592, 581, 569, 580],
[ 602, 592, 581, 572, 563, 579],
[ 579, 571, 563, 563, 564, 583],
[ 552, 551, 549, 562, 575, 592],
[ 537, 549, 561, 575, 590, 600]], dtype=int16)
amin[:6, :6] # minimum using a 3x3 moving window
array([[606, 592, 581, 569, 569, 569],
[592, 581, 572, 563, 563, 563],
[571, 563, 563, 563, 563, 563],
[551, 549, 549, 549, 562, 563],
[537, 537, 549, 549, 562, 564],
[537, 537, 549, 549, 562, 575]], dtype=uint16)
Now for the top 3x3 window, the minimum is 581 as can be seen by what is returned by the focal min. Do you want the location of the minimum? Since there are 2 locations for that one core cell as you will notice. Do you want the one on the top row, 3rd column or the 581 on the 3rd row, 3rd column... or both. It makes a difference, but is doable, if you can explain why that precision in location is needed.
So... if none of that applies, what is the number you want in the example above.
This thread has triggered my attention and I wanted to try something, just to get an idea of what we are talking about.
I created a python script that evaluates the largest drop in a moving window of 5x5 (still pretty small) to see what is possible, using a small extent of the data:
I converted the DEM into a point file and added the row and column numbers to it using the ($$RowMap and $$ColMap) in Python followed by the Extract Multivalues to Points tool.
Based on that data I created the script below where it travers through each center pixel (point) and extracts the values of the 5x5 neighborhood. It determines the min and max values and corresponding locations and generates a 3D polyline with attributes on drop (max - min) and some other stats.
If I color the lines by drop value and apply a minimum drop for the polyline to be generated (25m) this is the result:
I you look at a part of the area closer (showing only the largest drop values) you can see the effect:
Comparing it to a simple slope map:
and the lines on top of the slope raster
The code I used for this is listed below:
#-------------------------------------------------------------------------------
# Name: hydropower2.py
# Purpose:
#
# Author: Xander
#
# Created: 05-06-2017
#-------------------------------------------------------------------------------
def main():
import arcpy
import os
arcpy.env.overwriteOutput = True
fc = r'C:\GeoNet\HydroPower\data.gdb\points_dem_small'
fld_dem = 'grid_code'
fld_row = 'row'
fld_col = 'col'
fc_out = r'C:\GeoNet\HydroPower\data.gdb\drop_lines3'
fld_minh = "MinHeight"
fld_maxh = "MaxHeight"
fld_drop = "DropHeight"
fld_slope = "Slope"
# settings
min_drop = 25
# create a featureclass and fill with drop line
ws, fc_name = os.path.split(fc_out)
sr = arcpy.Describe(fc).spatialReference
arcpy.CreateFeatureclass_management(ws, fc_name, "POLYLINE", None, "DISABLED", "ENABLED", sr)
arcpy.AddField_management(fc_out, fld_minh, "DOUBLE")
arcpy.AddField_management(fc_out, fld_maxh, "DOUBLE")
arcpy.AddField_management(fc_out, fld_drop, "DOUBLE")
arcpy.AddField_management(fc_out, fld_slope, "DOUBLE")
flds = (fld_row, fld_col, fld_dem, 'SHAPE@')
dct = {'r{0}c{1}'.format(r[0], r[1]): (r[2], r[3]) for r in arcpy.da.SearchCursor(fc, flds)}
rows = max([r[0] for r in arcpy.da.SearchCursor(fc, (fld_row))])
cols = max([r[0] for r in arcpy.da.SearchCursor(fc, (fld_col))])
print "rows:", rows
print "cols:", cols
flds = ('SHAPE@', fld_minh, fld_maxh, fld_drop, fld_slope)
with arcpy.da.InsertCursor(fc_out, flds) as curs:
cnt = 0
drop_max = 0
for r in range(1, rows + 1):
for c in range(1, cols + 1):
cnt += 1
if cnt % 5000 == 0:
print "Processing", cnt
lst = getListNeighbors(dct, r, c, rows, cols)
min_val = getmin(lst)
max_val = getmax(lst)
drop = max_val - min_val
lst_pos_min = getPositions(lst, min_val)
lst_pos_max = getPositions(lst, max_val)
if drop >= min_drop:
min_d = 25
for pos_min in lst_pos_min:
for pos_max in lst_pos_max:
d = RelativeDistance(pos_min, pos_max)
if d < min_d:
min_d = d
pos_h = pos_max
pos_l = pos_min
polyline = CreatePolyline(lst, pos_h, pos_l)
curs.insertRow((polyline, min_val, max_val, drop, float(drop)/float(d), ))
if drop > drop_max:
drop_max = drop
print
print "drop_max:", drop_max
def CreatePolyline(lst, pos_h, pos_l):
tpl_h = GetItemOnTuple(lst, pos_h)
tpl_l = GetItemOnTuple(lst, pos_l)
pntg1 = tpl_h[1]
pntg2 = tpl_l[1]
sr = pntg1.spatialReference
pnt1 = pntg1.firstPoint
pnt1.Z = tpl_h[0]
pnt2 = pntg2.firstPoint
pnt2.Z = tpl_l[0]
return arcpy.Polyline(arcpy.Array([pnt1, pnt2]), sr, True, False)
def GetItemOnTuple(lst, tpl):
x = tpl[0]
y = tpl[1]
index = (-1*y+2)*5 + x + 2
return lst[index]
def RelativeDistance(tpl1, tpl2):
import math
return math.hypot(tpl1[0] - tpl2[0], tpl1[1] - tpl2[1])
def getmin(lst):
lst2 = [a[0] for a in lst if a[0] != None]
return min(lst2)
def getmax(lst):
lst2 = [a[0] for a in lst if a[0] != None]
return max(lst2)
def getPositions(lst, val):
positions = []
for r in range(5):
for c in range(5):
positions.append((c-2, -1*(r-2)))
result = []
for i in range(25):
if val == lst[i][0]:
result.append(positions[i])
return result
def getListNeighbors(dct, ri, ci, rows, cols):
lst = []
prn = True
for r in range(ri - 2, ri + 3):
for c in range(ci - 2, ci + 3):
rc = 'r{0}c{1}'.format(r, c)
if rc in dct:
lst.append(dct[rc])
else:
prn = False
lst.append((None, None))
return lst
if __name__ == '__main__':
main()You could apply this to larger moving windows, but I'm still not convinced by it...
