In [ ]:
%pylab inline
from PIL import Image
from scipy.signal import convolve2d
# Neuer tif-Link https://my.mail.de/dl/8c55544751dd4272a1c558a98fa545e2
Disclaimer: alles folgende wurde nach bestem Gewissen erstellt. Fehler dürfen gefunden und korrigiert werden.
In [ ]:
terrain = np.array(Image.open("srtm_39_02b.tif"))
In [ ]:
# Table 1, Step 1
# Wir füllen einzelne Mulden so hoch wie die niedrigste umgebende Zelle hoch ist
def fill_depressions(original_terrain):
terrain = np.copy(original_terrain)
for i in xrange(1, terrain.shape[0]-1):
if i % 100 == 0:
print i
for j in xrange(1, terrain.shape[1]-1):
# Ränder werden erstmal weggelassen
# wir betrachten Pixel (i, j)
pixel_is_depr = True
min_height = 10e5 # irgendein sehr hoher Wert
for k in xrange(-1, 2):
for l in range(-1, 2):
if (i, j) != (0, 0):
# sobald ein nachbar nicht höher als pixel i,j ist, ist es keine Senke
pixel_is_depr = pixel_is_depr and terrain[i+k, j+l] > terrain[i, j]
min_height = min(min_height, terrain[i+k, j+l])
if pixel_is_depr:
terrain[i, j] = min_height
return terrain
# sehr langsam!
In [ ]:
# Numpy-optimierte Version
def fill_dep_better(original_terrain):
#terrain = np.copy(original_terrain)
# Wir bauen einen Array, der aus 9 "Schichten" besteht:
# Erste Schicht: Originalarray ohne Ränder
# Jede weitere Schicht um eins horizontal, vertikal oder diagonal verschoben
schicht_array = np.zeros((9, original_terrain.shape[0]-2, original_terrain.shape[1]-2), int)
schicht_array[0] = original_terrain[1:-1, 1:-1]
schicht = 1
for k in xrange(-1, 2):
for l in xrange(-1, 2):
if (k, l) != (0, 0):
#print np.roll(np.roll(original_terrain, k, axis=0), l, axis=1)
schicht_array[schicht] = np.roll(np.roll(original_terrain, k, axis=0), l, axis=1)[1:-1, 1:-1]
schicht += 1
# entlang Achse 0 liegt jetzt jeder Pixel in schicht_array[0] mit seinen Nachbarn übereinander
res = (schicht_array[1:] > schicht_array[0]).all(axis=0)
ind1, ind2 = np.where(res)
schicht_array[0, ind1, ind2] = schicht_array[1:, ind1, ind2].min(axis=0)
return schicht_array[0], res
#Image.fromarray(np.array(res, dtype=np.uint8)).show()
In [ ]:
t2, res = fill_dep_better(terrain)
In [ ]:
indices = zip(*np.where(t2 != terrain[1:-1, 1:-1]))
In [ ]:
for i, j in indices:
print terrain[i:i+3, j:j+3]
print t2[i-1:i+2, j-1:j+2]
In [ ]:
terrain.shape
In [ ]: