

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.

Schritt 1 : Mulden (Depressions) füllen

Wasserfluss lässt sich schwierig von DEMs bestimmen, wenn Mulden vorhanden sind -> Mulden füllen

# Flussrichtungen bestimmen # Jede Zelle hat 8 mögliche Flussrichtungen Diese sind folgendermaßen kodiert: 64 128 1 32 x 2 16 8 4



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 [ ]: