csmap-pyを動く状態でQGISプラグインに統合

__init__.py

@@ -3,9 +3,10 @@
 
 # to import modules as non-relative
 sys.path.append(os.path.dirname(__file__))
+sys.path.append(os.path.join(os.path.dirname(__file__), "csmap_py"))
 
 
 def classFactory(iface):
-    from csmap import CsMap
+    from qcsmap import QCsMap
 
-    return CsMap(iface)
+    return QCsMap(iface)

csmap_py/csmap/__main__.py

@@ -0,0 +1,82 @@
+from csmap.process import CsmapParams, process
+
+
+def parse_args():
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("input_dem_path", type=str, help="input DEM path")
+    parser.add_argument("output_path", type=str, help="output path")
+    parser.add_argument(
+        "--chunk_size",
+        type=int,
+        default=1024,
+        help="chunk size as pixel, default to 1024",
+    )
+    parser.add_argument(
+        "--max_workers",
+        type=int,
+        default=1,
+        help="max workers for multiprocessing, default to 1",
+    )
+    parser.add_argument(
+        "--gf_size", type=int, default=12, help="gaussian filter size, default to 12"
+    )
+    parser.add_argument(
+        "--gf_sigma", type=int, default=3, help="gaussian filter sigma, default to 3"
+    )
+    parser.add_argument(
+        "--curvature_size",
+        type=int,
+        default=1,
+        help="curvature filter size, default to 1",
+    )
+    parser.add_argument(
+        "--height_scale",
+        type=float,
+        nargs=2,
+        default=[0.0, 1000.0],
+        help="height scale, min max, default to 0.0 1000.0",
+    )
+    parser.add_argument(
+        "--slope_scale",
+        type=float,
+        nargs=2,
+        default=[0.0, 1.5],
+        help="slope scale, min max, default to 0.0 1.5",
+    )
+    parser.add_argument(
+        "--curvature_scale",
+        type=float,
+        nargs=2,
+        default=[-0.1, 0.1],
+        help="curvature scale, min max, default to -0.1 0.1",
+    )
+
+    args = parser.parse_args()
+
+    params = CsmapParams(
+        gf_size=args.gf_size,
+        gf_sigma=args.gf_sigma,
+        curvature_size=args.curvature_size,
+        height_scale=args.height_scale,
+        slope_scale=args.slope_scale,
+        curvature_scale=args.curvature_scale,
+    )
+
+    return {
+        "input_dem_path": args.input_dem_path,
+        "output_path": args.output_path,
+        "chunk_size": args.chunk_size,
+        "max_workers": args.max_workers,
+        "params": params,
+    }
+
+
+def main():
+    args = parse_args()
+    process(**args)
+
+
+if __name__ == "__main__":
+    main()

csmap_py/csmap/calc.py

@@ -0,0 +1,75 @@
+import numpy as np
+
+
+def slope(dem: np.ndarray) -> np.ndarray:
+    """傾斜率を求める
+    出力のndarrayのshapeは、(dem.shape[0] - 2, dem.shape[1] - 2)
+    """
+    z2 = dem[1:-1, 0:-2]
+    z4 = dem[0:-2, 1:-1]
+    z6 = dem[2:, 1:-1]
+    z8 = dem[1:-1, 2:]
+    p = (z6 - z4) / 2
+    q = (z8 - z2) / 2
+    p2 = p * p
+    q2 = q * q
+
+    slope = np.arctan((p2 + q2) ** 0.5)
+    return slope
+
+
+def gaussianfilter(image: np.ndarray, size: int, sigma: int) -> np.ndarray:
+    """ガウシアンフィルター"""
+    size = int(size) // 2
+    x, y = np.mgrid[-size : size + 1, -size : size + 1]
+    g = np.exp(-(x**2 + y**2) / (2 * sigma**2))
+    kernel = g / g.sum()
+
+    # 画像を畳み込む
+    k_h, k_w = kernel.shape
+    i_h, i_w = image.shape
+
+    # パディングサイズを計算
+    pad_h, pad_w = k_h // 2, k_w // 2
+
+    # 画像にパディングを適用
+    padded = np.pad(image, ((pad_h, pad_h), (pad_w, pad_w)), mode="constant")
+
+    # einsumを使用して畳み込みを行う
+    sub_matrices = np.lib.stride_tricks.as_strided(
+        padded, shape=(i_h, i_w, k_h, k_w), strides=padded.strides * 2
+    )
+    return np.einsum("ijkl,kl->ij", sub_matrices, kernel)
+
+
+def curvature(dem: np.ndarray, cell_size: int) -> np.ndarray:
+    """曲率を求める"""
+
+    # SAGA の Slope, Aspect, Curvature の 9 parameter 2nd order polynom に準拠
+    z1 = dem[0:-2, 0:-2]
+    z2 = dem[1:-1, 0:-2]
+    z3 = dem[2:, 0:-2]
+    z4 = dem[0:-2, 1:-1]
+    z5 = dem[1:-1, 1:-1]
+    z6 = dem[2:, 1:-1]
+    z7 = dem[0:-2, 2:]
+    z8 = dem[1:-1, 2:]
+    z9 = dem[2:, 2:]
+
+    cell_area = cell_size * cell_size
+    r = ((z4 + z6) / 2 - z5) / cell_area
+    t = ((z2 + z8) / 2 - z5) / cell_area
+    p = (z6 - z4) / (2 * cell_size)
+    q = (z8 - z2) / (2 * cell_size)
+    s = (z1 - z3 - z7 + z9) / (4 * cell_area)
+    p2 = p * p
+    q2 = q * q
+    spq = s * p * q
+
+    # gene
+    return -2 * (r + t)
+
+    # plan
+    p2_q2 = p2 + q2
+    p2_q2 = np.where(p2_q2 > 1e-6, p2_q2, np.nan)
+    return -(t * p2 + r * q2 - 2 * spq) / ((p2_q2) ** 1.5)

csmap_py/csmap/color.py

@@ -0,0 +1,112 @@
+import numpy as np
+
+
+def rgbify(arr: np.ndarray, method, scale: (float, float) = None) -> np.ndarray:
+    """ndarrayをRGBに変換する
+    - arrは変更しない
+    - ndarrayのshapeは、(4, height, width) 4はRGBA
+    """
+
+    _min = arr.min() if scale is None else scale[0]
+    _max = arr.max() if scale is None else scale[1]
+
+    # -x ~ x を 0 ~ 1 に正規化
+    arr = (arr - _min) / (_max - _min)
+    # clamp
+    arr = np.where(arr < 0, 0, arr)
+    arr = np.where(arr > 1, 1, arr)
+
+    # 3次元に変換
+    rgb = method(arr)
+    return rgb
+
+
+def slope_red(arr: np.ndarray) -> np.ndarray:
+    rgb = np.zeros((4, arr.shape[0], arr.shape[1]), dtype=np.uint8)
+    rgb[0, :, :] = 255 - arr * 155  # R: 255 -> 100
+    rgb[1, :, :] = 245 - arr * 195  # G: 245 -> 50
+    rgb[2, :, :] = 235 - arr * 215  # B: 235 -> 20
+    rgb[3, :, :] = 255
+    return rgb
+
+
+def slope_blackwhite(arr: np.ndarray) -> np.ndarray:
+    rgb = np.zeros((4, arr.shape[0], arr.shape[1]), dtype=np.uint8)
+    rgb[0, :, :] = (1 - arr) * 255  # R
+    rgb[1, :, :] = (1 - arr) * 255  # G
+    rgb[2, :, :] = (1 - arr) * 255  # B
+    rgb[3, :, :] = 255  # A
+    return rgb
+
+
+def curvature_blue(arr: np.ndarray) -> np.ndarray:
+    rgb = np.zeros((4, arr.shape[0], arr.shape[1]), dtype=np.uint8)
+    rgb[0, :, :] = 35 + arr * 190  # R: 35 -> 225
+    rgb[1, :, :] = 80 + arr * 155  # G: 80 -> 235
+    rgb[2, :, :] = 100 + arr * 145  # B: 100 -> 245
+    rgb[3, :, :] = 255
+    return rgb
+
+
+def curvature_redyellowblue(arr: np.ndarray) -> np.ndarray:
+    # value:0-1 to: red -> yellow -> blue
+    # interpolate between red and yellow, and yellow and blue, by linear
+
+    # 0-0.5: blue -> white
+    rgb1 = np.zeros((4, arr.shape[0], arr.shape[1]), dtype=np.uint8)
+    rgb1[0, :, :] = 75 + arr * 170 * 2  # R: 75 -> 245
+    rgb1[1, :, :] = 100 + arr * 145 * 2  # G: 100 -> 245
+    rgb1[2, :, :] = 165 + arr * 80 * 2  # B: 165 -> 245
+
+    # 0.5-1: white -> red
+    rgb2 = np.zeros((4, arr.shape[0], arr.shape[1]), dtype=np.uint8)
+    rgb2[0, :, :] = 245 - (arr * 2 - 1) * 100  # R: 245 -> 145
+    rgb2[1, :, :] = 245 - (arr * 2 - 1) * 190  # G: 245 -> 55
+    rgb2[2, :, :] = 245 - (arr * 2 - 1) * 195  # B: 245 -> 50
+
+    # blend
+    rgb = np.where(arr < 0.5, rgb1, rgb2)
+    rgb[3, :, :] = 255
+
+    return rgb
+
+
+def height_blackwhite(arr: np.ndarray) -> np.ndarray:
+    rgb = np.zeros((4, arr.shape[0], arr.shape[1]), dtype=np.uint8)
+    rgb[0, :, :] = (1 - arr) * 255  # R
+    rgb[1, :, :] = (1 - arr) * 255  # G
+    rgb[2, :, :] = (1 - arr) * 255  # B
+    rgb[3, :, :] = 255
+    return rgb
+
+
+def blend(
+    dem_bw: np.ndarray,
+    slope_red: np.ndarray,
+    slope_bw: np.ndarray,
+    curvature_blue: np.ndarray,
+    curvature_ryb: np.ndarray,
+    blend_params: dict = {
+        "slope_bw": 0.5,  # alpha blending based on the paper
+        "curvature_ryb": 0.25,  # 0.5 / 2
+        "slope_red": 0.125,  # 0.5 / 2 / 2
+        "curvature_blue": 0.06125,  # 0.5 / 2 / 2 / 2
+        "dem": 0.030625,  # 0.5 / 2 / 2 / 2 / 2
+    },
+) -> np.ndarray:
+    """blend all rgb
+    全てのndarrayは同じshapeであること
+    DEMを用いて処理した他の要素は、DEMよりも1px内側にpaddingされているので
+    あらかじめDEMのpaddingを除外しておく必要がある
+    """
+    _blend = np.zeros((4, dem_bw.shape[0], dem_bw.shape[1]), dtype=np.uint8)
+    _blend = (
+        dem_bw * blend_params["dem"]
+        + slope_red * blend_params["slope_red"]
+        + slope_bw * blend_params["slope_bw"]
+        + curvature_blue * blend_params["curvature_blue"]
+        + curvature_ryb * blend_params["curvature_ryb"]
+    )
+    _blend = _blend.astype(np.uint8)  # force uint8
+    _blend[3, :, :] = 255  # alpha
+    return _blend