From 884d6cd57d4aa8e79ba560e17959088c42f72d35 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 20:59:23 +0900 Subject: [PATCH 1/9] =?UTF-8?q?fix:=20csmap=5Fpy=E3=82=92PYTHONPATH?= =?UTF-8?q?=E3=81=AB=E8=BF=BD=E5=8A=A0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- __init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/__init__.py b/__init__.py index 1aef0e6..deb4cd2 100644 --- a/__init__.py +++ b/__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) From 17313185da990c4c6d5541cf8eee9f913d79b1e5 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 20:59:48 +0900 Subject: [PATCH 2/9] =?UTF-8?q?fix:=20csmap-py=E3=81=A8=E5=90=8D=E5=89=8D?= =?UTF-8?q?=E3=81=8C=E8=A1=9D=E7=AA=81=E3=81=97=E3=81=A6=E3=81=84=E3=82=8B?= =?UTF-8?q?=E3=81=AE=E3=81=A7qcsmap=E3=81=AB=E3=83=AA=E3=83=8D=E3=83=BC?= =?UTF-8?q?=E3=83=A0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- qcsmap.py | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 qcsmap.py diff --git a/qcsmap.py b/qcsmap.py new file mode 100644 index 0000000..5ced4a6 --- /dev/null +++ b/qcsmap.py @@ -0,0 +1,65 @@ +import os + +from PyQt5.QtWidgets import QAction +from qgis.PyQt.QtGui import QIcon + +from dem_to_csmap import DemToCsMap + +PLUGIN_NAME = "CSMap Plugin" + + +class QCsMap: + def __init__(self, iface): + self.iface = iface + self.win = self.iface.mainWindow() + self.plugin_dir = os.path.dirname(__file__) + self.actions = [] + self.menu = PLUGIN_NAME + self.toolbar = self.iface.addToolBar(PLUGIN_NAME) + self.toolbar.setObjectName(PLUGIN_NAME) + + def add_action( + self, + icon_path, + text, + callback, + enabled_flag=True, + add_to_menu=True, + add_to_toolbar=True, + status_tip=None, + whats_this=None, + parent=None, + ): + icon = QIcon(icon_path) + action = QAction(icon, text, parent) + action.triggered.connect(callback) + action.setEnabled(enabled_flag) + if status_tip is not None: + action.setStatusTip(status_tip) + if whats_this is not None: + action.setWhatsThis(whats_this) + if add_to_toolbar: + self.toolbar.addAction(action) + if add_to_menu: + self.iface.addPluginToMenu(self.menu, action) + self.actions.append(action) + return action + + def initGui(self): + # アイコンのパスを取得 + icon_path = os.path.join(os.path.dirname(__file__), "imgs", "icon.png") + # メニュー設定 + self.add_action( + icon_path=icon_path, text="", callback=self.show_menu_dem_to_csmap, parent=self.win + ) + + def unload(self): + for action in self.actions: + self.iface.removePluginMenu(PLUGIN_NAME, action) + self.iface.removeToolBarIcon(action) + del self.toolbar + + def show_menu_dem_to_csmap(self): + self.show_menu_dem_to_csmap = DemToCsMap() + self.show_menu_dem_to_csmap.show() + From c65aa2e1f1a848532507cf4966b651637ca9abc8 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:00:10 +0900 Subject: [PATCH 3/9] =?UTF-8?q?sample:=20csmap-py=E3=82=92=E5=AE=9F?= =?UTF-8?q?=E8=A1=8C=E3=81=99=E3=82=8B=E3=82=B5=E3=83=B3=E3=83=97=E3=83=AB?= =?UTF-8?q?=E3=82=B3=E3=83=BC=E3=83=89?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- dem_to_csmap.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/dem_to_csmap.py b/dem_to_csmap.py index 890e9a5..fda00e4 100644 --- a/dem_to_csmap.py +++ b/dem_to_csmap.py @@ -3,6 +3,8 @@ from PyQt5.QtWidgets import QDialog, QMessageBox from qgis.PyQt import uic +from csmap_py.csmap import process + class DemToCsMap(QDialog): def __init__(self): @@ -19,3 +21,14 @@ def get_and_show_input_text(self): text_value = self.ui.lineEdit.text() # テキストボックス値をメッセージ表示 QMessageBox.information(None, "ウィンドウ名", text_value) + + # 試しにCSMapの処理を実行：ちゃんと入力・出力をUIから参照しよう + params = process.CsmapParams() + input_path = "/Users/kanahiro/Downloads/dem.tif" + output_path = "/Users/kanahiro/Downloads/out.tif" + process.process( + input_path, + output_path, + 256, + params=params, + ) From 0880be0fe2105746696213aef96fbf883cc60628 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:01:08 +0900 Subject: [PATCH 4/9] feat: add csmap-py --- csmap_py | 1 + 1 file changed, 1 insertion(+) create mode 160000 csmap_py diff --git a/csmap_py b/csmap_py new file mode 160000 index 0000000..a9bf18a --- /dev/null +++ b/csmap_py @@ -0,0 +1 @@ +Subproject commit a9bf18a6c9e354707d4c410dd2b0d7fdbd35b01b From 80571b5f84e34e9fa51c0c4533cf6e757dd909e1 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:09:26 +0900 Subject: [PATCH 5/9] =?UTF-8?q?=E4=B8=80=E6=97=A6csmap-py=E6=B6=88?= =?UTF-8?q?=E3=81=99?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- csmap.py | 65 -------------------------------------------------------- csmap_py | 1 - 2 files changed, 66 deletions(-) delete mode 100644 csmap.py delete mode 160000 csmap_py diff --git a/csmap.py b/csmap.py deleted file mode 100644 index f24b52a..0000000 --- a/csmap.py +++ /dev/null @@ -1,65 +0,0 @@ -import os - -from PyQt5.QtWidgets import QAction -from qgis.PyQt.QtGui import QIcon - -from dem_to_csmap import DemToCsMap - -PLUGIN_NAME = "CSMap Plugin" - - -class CsMap: - def __init__(self, iface): - self.iface = iface - self.win = self.iface.mainWindow() - self.plugin_dir = os.path.dirname(__file__) - self.actions = [] - self.menu = PLUGIN_NAME - self.toolbar = self.iface.addToolBar(PLUGIN_NAME) - self.toolbar.setObjectName(PLUGIN_NAME) - - def add_action( - self, - icon_path, - text, - callback, - enabled_flag=True, - add_to_menu=True, - add_to_toolbar=True, - status_tip=None, - whats_this=None, - parent=None, - ): - icon = QIcon(icon_path) - action = QAction(icon, text, parent) - action.triggered.connect(callback) - action.setEnabled(enabled_flag) - if status_tip is not None: - action.setStatusTip(status_tip) - if whats_this is not None: - action.setWhatsThis(whats_this) - if add_to_toolbar: - self.toolbar.addAction(action) - if add_to_menu: - self.iface.addPluginToMenu(self.menu, action) - self.actions.append(action) - return action - - def initGui(self): - # アイコンのパスを取得 - icon_path = os.path.join(os.path.dirname(__file__), "imgs", "icon.png") - # メニュー設定 - self.add_action( - icon_path=icon_path, text="", callback=self.show_menu_dem_to_csmap, parent=self.win - ) - - def unload(self): - for action in self.actions: - self.iface.removePluginMenu(PLUGIN_NAME, action) - self.iface.removeToolBarIcon(action) - del self.toolbar - - def show_menu_dem_to_csmap(self): - self.show_menu_dem_to_csmap = DemToCsMap() - self.show_menu_dem_to_csmap.show() - diff --git a/csmap_py b/csmap_py deleted file mode 160000 index a9bf18a..0000000 --- a/csmap_py +++ /dev/null @@ -1 +0,0 @@ -Subproject commit a9bf18a6c9e354707d4c410dd2b0d7fdbd35b01b From ba049e07079aaabb3b1823b496daf69e9c45e7ce Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:10:08 +0900 Subject: [PATCH 6/9] add csmap-py --- csmap_py/csmap/__init__.py | 0 csmap_py/csmap/__main__.py | 82 +++++++++++++++++ csmap_py/csmap/calc.py | 75 +++++++++++++++ csmap_py/csmap/color.py | 112 +++++++++++++++++++++++ csmap_py/csmap/process.py | 182 +++++++++++++++++++++++++++++++++++++ 5 files changed, 451 insertions(+) create mode 100644 csmap_py/csmap/__init__.py create mode 100644 csmap_py/csmap/__main__.py create mode 100644 csmap_py/csmap/calc.py create mode 100644 csmap_py/csmap/color.py create mode 100644 csmap_py/csmap/process.py diff --git a/csmap_py/csmap/__init__.py b/csmap_py/csmap/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/csmap_py/csmap/__main__.py b/csmap_py/csmap/__main__.py new file mode 100644 index 0000000..9a4d62b --- /dev/null +++ b/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() diff --git a/csmap_py/csmap/calc.py b/csmap_py/csmap/calc.py new file mode 100644 index 0000000..485509f --- /dev/null +++ b/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) diff --git a/csmap_py/csmap/color.py b/csmap_py/csmap/color.py new file mode 100644 index 0000000..c8f67d2 --- /dev/null +++ b/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 diff --git a/csmap_py/csmap/process.py b/csmap_py/csmap/process.py new file mode 100644 index 0000000..e77f2be --- /dev/null +++ b/csmap_py/csmap/process.py @@ -0,0 +1,182 @@ +from dataclasses import dataclass +from threading import Lock +from concurrent import futures + +import numpy as np +import rasterio +from rasterio.windows import Window +from rasterio.transform import Affine + +from csmap import calc +from csmap import color + + +@dataclass +class CsmapParams: + gf_size: int = 12 + gf_sigma: int = 3 + curvature_size: int = 1 + height_scale: (float, float) = (0.0, 1000.0) + slope_scale: (float, float) = (0.0, 1.5) + curvature_scale: (float, float) = (-0.1, 0.1) + + +def csmap(dem: np.ndarray, params: CsmapParams) -> np.ndarray: + """DEMからCS立体図を作成する""" + # calclucate elements + slope = calc.slope(dem) + g = calc.gaussianfilter(dem, params.gf_size, params.gf_sigma) + curvature = calc.curvature(g, params.curvature_size) + + # rgbify + dem_rgb = color.rgbify(dem, color.height_blackwhite, scale=params.height_scale) + slope_red = color.rgbify(slope, color.slope_red, scale=params.slope_scale) + slope_bw = color.rgbify(slope, color.slope_blackwhite, scale=params.slope_scale) + curvature_blue = color.rgbify( + curvature, color.curvature_blue, scale=params.curvature_scale + ) + curvature_ryb = color.rgbify( + curvature, color.curvature_redyellowblue, scale=params.curvature_scale + ) + + dem_rgb = dem_rgb[:, 1:-1, 1:-1] # remove padding + + # blend all rgb + blend = color.blend( + dem_rgb, + slope_red, + slope_bw, + curvature_blue, + curvature_ryb, + ) + + return blend + + +def _process_chunk( + chunk: np.ndarray, + dst: rasterio.io.DatasetWriter, + x: int, + y: int, + write_size_x: int, + write_size_y: int, + params: CsmapParams, + lock: Lock = None, +) -> np.ndarray: + """チャンクごとの処理""" + csmap_chunk = csmap(chunk, params) + csmap_chunk_margin_removed = csmap_chunk[ + :, + (params.gf_size + params.gf_sigma) + // 2 : -((params.gf_size + params.gf_sigma) // 2), + (params.gf_size + params.gf_sigma) + // 2 : -((params.gf_size + params.gf_sigma) // 2), + ] # shape = (4, chunk_size - margin, chunk_size - margin) + + if lock is None: + dst.write( + csmap_chunk_margin_removed, + window=Window(x, y, write_size_x, write_size_y), + ) + else: + with lock: + dst.write( + csmap_chunk_margin_removed, + window=Window(x, y, write_size_x, write_size_y), + ) + + +def process( + input_dem_path: str, + output_path: str, + chunk_size: int, + params: CsmapParams, + max_workers: int = 1, +): + with rasterio.open(input_dem_path) as dem: + margin = params.gf_size + params.gf_sigma # ガウシアンフィルタのサイズ+シグマ + # チャンクごとの処理結果には「淵=margin」が生じるのでこの部分を除外する必要がある + margin_to_removed = margin // 2 # 整数値に切り捨てた値*両端 + + # マージンを考慮したtransform + transform = Affine( + dem.transform.a, + dem.transform.b, + dem.transform.c + (1 + margin // 2) * dem.transform.a, # 左端の座標をマージン分ずらす + dem.transform.d, + dem.transform.e, + dem.transform.f + (1 + margin // 2) * dem.transform.e, # 上端の座標をマージン分ずらす + 0.0, + 0.0, + 1.0, + ) + + # 生成されるCS立体図のサイズ + out_width = dem.shape[1] - margin_to_removed * 2 - 2 + out_height = dem.shape[0] - margin_to_removed * 2 - 2 + + with rasterio.open( + output_path, + "w", + driver="GTiff", + dtype=rasterio.uint8, + count=4, + width=out_width, + height=out_height, + crs=dem.crs, + transform=transform, + compress="LZW", + ) as dst: + # chunkごとに処理 + chunk_csmap_size = chunk_size - margin_to_removed * 2 - 2 + + # 並列処理しない場合とする場合で処理を分ける + if max_workers == 1: + for y in range(0, dem.shape[0], chunk_csmap_size): + for x in range(0, dem.shape[1], chunk_csmap_size): + # csmpのどの部分を出力用配列に入れるかを計算 + write_size_x = chunk_csmap_size + write_size_y = chunk_csmap_size + if x + chunk_csmap_size > out_width: + write_size_x = out_width - x + if y + chunk_csmap_size > out_height: + write_size_y = out_height - y + + chunk = dem.read(1, window=Window(x, y, chunk_size, chunk_size)) + _process_chunk( + chunk, + dst, + x, + y, + write_size_x, + write_size_y, + params, + ) + else: # 並列処理する場合=ThreadPoolExecutorを使用する + lock = Lock() # 並列処理のロック + with futures.ThreadPoolExecutor(max_workers=max_workers) as executor: + # chunkごとに処理 + for y in range(0, dem.shape[0], chunk_csmap_size): + for x in range(0, dem.shape[1], chunk_csmap_size): + # csmpのどの部分を出力用配列に入れるかを計算 + write_size_x = chunk_csmap_size + write_size_y = chunk_csmap_size + if x + chunk_csmap_size > out_width: + write_size_x = out_width - x + if y + chunk_csmap_size > out_height: + write_size_y = out_height - y + + chunk = dem.read( + 1, window=Window(x, y, chunk_size, chunk_size) + ) + executor.submit( + _process_chunk, + chunk, + dst, + x, + y, + write_size_x, + write_size_y, + params, + lock, + ) From 0877ae5e2d7dfa190d50255a9e8b9e4dae68d915 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:10:53 +0900 Subject: [PATCH 7/9] add csmap-py --- csmap_py/csmap/process.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/csmap_py/csmap/process.py b/csmap_py/csmap/process.py index e77f2be..ed029a6 100644 --- a/csmap_py/csmap/process.py +++ b/csmap_py/csmap/process.py @@ -106,10 +106,7 @@ def process( dem.transform.d, dem.transform.e, dem.transform.f + (1 + margin // 2) * dem.transform.e, # 上端の座標をマージン分ずらす - 0.0, - 0.0, - 1.0, - ) + ) # 古いrasterioに合わせて引数を3つ削除した # 生成されるCS立体図のサイズ out_width = dem.shape[1] - margin_to_removed * 2 - 2 From d3a32714d64bbec7f59d32a381eee6d22c04d15c Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:18:54 +0900 Subject: [PATCH 8/9] add rasterio in dev --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index df0b921..147b070 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,7 @@ python = "^3.9" [tool.poetry.group.dev.dependencies] ruff = "^0.4.5" +rasterio = "^1.3.10" [build-system] From b0d7352f14a7131fff532fe7f8dda6ceb400d1a3 Mon Sep 17 00:00:00 2001 From: Kanahiro Date: Wed, 5 Jun 2024 21:19:09 +0900 Subject: [PATCH 9/9] =?UTF-8?q?fix:=20csmap=E3=82=92import=E3=81=99?= =?UTF-8?q?=E3=82=8B=E3=81=9F=E3=82=81=E3=81=AE=E5=87=A6=E7=BD=AE?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- tests/__init__.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/__init__.py b/tests/__init__.py index e69de29..63d8e0e 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -0,0 +1,5 @@ +import os +import sys + +# csmap_pyからcsmapとしてimportするためにパスを追加 +sys.path.append(os.path.join(os.path.dirname(os.path.dirname(__file__)), "csmap_py"))