TopoAnalysis/demMethods.py at master · stgl/TopoAnalysis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
"""Batch pipelines that chain the grid classes into a finished analysis."""

import numpy as np

try:  # pragma: no cover - exercised by whichever import style is in use
    from . import dem as d
except ImportError:  # pragma: no cover
    import dem as d

def processAll(prefix_name, Ao, theta, base_name = '.'):

    elevation_name = base_name + "/" + prefix_name + "_dem_15s"
    area_name = base_name + "/" + prefix_name + "_acc_15s"
    d8_name = base_name + "/" + prefix_name + "_dir_15s"

    elevation = d.Elevation(gdal_filename = elevation_name)
    elevation.save(prefix_name + "_elevation")
    area = d.Area(gdal_filename = area_name)
    d8 = d.FlowDirectionD8(gdal_filename = d8_name)
    d8.save(prefix_name + "_flow_direction")

    idx = area.sort(reverse = False)
    area = d.GeographicArea(flow_direction = d8, sorted_indexes = idx)
    area.save(prefix_name + "_area")
    flow_length = d.GeographicFlowLength(flow_direction = d8, sorted_indexes = idx)
    flow_length.save(prefix_name + "_flow_length")
    ksi = d.GeographicKsi(area = area, flow_direction = d8, theta = theta, Ao = Ao, flow_length = flow_length, sorted_indexes = idx)
    ksi.save(prefix_name + "_ksi_" + str(Ao).replace('.','_') + "_" + str(theta).replace('.','_'))
    relief = d.ScaledRelief(flow_direction = d8, elevation = elevation, flow_length = flow_length, Ao = Ao, theta = theta, sorted_indexes = idx)
    relief.save(prefix_name + "_relief_" + str(Ao).replace('.','_') + "_" + str(theta).replace('.','_'))

def processAllUTM(prefix_name, EPSGprojectionCode, Ao, theta, base_name = '.'):

    full_path_without_suffix = base_name + "/" + prefix_name
    elevation_unfilled_ascii_filename = full_path_without_suffix + ".txt"
    elevation = d.Elevation(ai_ascii_filename = elevation_unfilled_ascii_filename, EPSGprojectionCode= EPSGprojectionCode)
    elevation.save(full_path_without_suffix + "_elevation")
    filled = d.FilledElevation(elevation = elevation)
    filled.save(full_path_without_suffix + "_filled")
    d8 = d.FlowDirectionD8(flooded_dem = filled)
    d8.save(full_path_without_suffix + "_flow_direction")
    area = d.Area(flow_direction = d8)
    area.save(full_path_without_suffix + "_area")

    flow_length = d.FlowLength(flow_direction = d8)
    flow_length.save(full_path_without_suffix + "_flow_length")
    ksi = d.Ksi(area = area, flow_direction = d8, theta = theta, Ao = Ao, flow_length = flow_length)
    ksi.save(full_path_without_suffix + "_ksi_" + str(Ao).replace('.','_') + "_" + str(theta).replace('.','_'))
    relief = d.ScaledRelief(flow_direction = d8, elevation = elevation, flow_length = flow_length, Ao = Ao, theta = theta)
    relief.save(full_path_without_suffix + "_relief_" + str(Ao).replace('.','_') + "_" + str(theta).replace('.','_'))

def processForTheta(prefix_name, Ao, theta, base_name = '.'):

    area_name = base_name + "/" + prefix_name + "_area"
    d8_name = base_name + "/" + prefix_name + "_flow_direction"
    flow_length_name = base_name + "/" + prefix_name + "_flow_length"
    elevation_name = base_name + "/" + prefix_name + "_elevation"

    area = d.GeographicArea.load(area_name)
    d8 = d.FlowDirectionD8.load(d8_name)
    flow_length = d.GeographicFlowLength.load(flow_length_name)
    elevation = d.Elevation.load(elevation_name)
    idx = area.sort(reverse = False)
    ksi = d.GeographicKsi(area = area, flow_direction = d8, theta = theta, Ao = Ao, flow_length = flow_length, sorted_indexes = idx)
    ksi.save(prefix_name + "_ksi_" + str(Ao).replace('.','_') + "_" + str(theta).replace('.','_'))
    relief = d.ScaledRelief(flow_direction = d8, elevation = elevation, flow_length = flow_length, Ao = Ao, theta = theta, sorted_indexes = idx, area = area)
    relief.save(prefix_name + "_relief_" + str(Ao).replace('.','_') + "_" + str(theta).replace('.','_'))

def plotGrids(x_grid, y_grid, plot_string, **kwargs):
    """Scatter one grid against another, keeping only non-negative y."""
    x_vec = np.ndarray.flatten(x_grid._griddata)
    y_vec = np.ndarray.flatten(y_grid._griddata)

    ind = np.where(y_vec >= 0)
    x_vec = x_vec[ind]
    y_vec = y_vec[ind]

    from matplotlib import pyplot as plt

    plt.plot(x_vec, y_vec, plot_string, **kwargs)

    return x_vec, y_vec

def extract_values_from_grid(x_grid, y_grid, ignore_zeros=False):
    """Flatten two grids to matching 1-D vectors."""
    x_vec = np.ndarray.flatten(x_grid._griddata)
    y_vec = np.ndarray.flatten(y_grid._griddata)
    if ignore_zeros:
        i = np.where(np.logical_and(x_vec != 0, y_vec != 0))
        x_vec = x_vec[i[0]]
        y_vec = y_vec[i[0]]

    return x_vec, y_vec

def create_density(x, y, x_boundaries, y_boundaries):
    """2-D histogram of paired grid values."""
    H, xedges, yedges = np.histogram2d(x, y, bins = (x_boundaries, y_boundaries))

    return H, xedges, yedges