池の処理（１）
池の処理を考えてみる。
以下の処理は、3次メッシュを対象として考える。
用意するデータ
- 3次メッシュの標高データ（data_array）
- 3次メッシュのDEM構成点種別マップ（surf_array）
- 上のDEM構成点種別マップから物体の数を数えるの結果のlabelデータ（obj_array）
これら3つのデータをロードした状態で以下の処理を行う。
- obj_arrayから対象とする物体（池）の番号を指定して、範囲を求める。
- 上で求めた範囲について、obj_array上で池がある場合に、その地点のsurf_arrayでの上下左右の地点が「地表面」の場所の標高をdata_arrayから取得し、最低標高を求める。（池の岸の標高の最低地点を池の水面の標高とする。）
- 池の底の標高を池の水面の標高-5mとして、obj_arrayの範囲内を対象にして、池の部分に池の底の標高をセットする。
- data_arrayを出力する。
以上の処理をobjlake.pyとして書いて見た。
# coding: utf-8 import numpy as np import re import os import sys NPY_DIR = "/Users/pycra/Desktop/NPY_DATA" SURF_FILENAME = "npysurf-{0:04d}-{1:02d}-{2:02d}.npy" DATA_FILENAME = "npydata-{0:04d}-{1:02d}-{2:02d}.npy" TEMP_LAKE_FILENAME = "temp_lake.npy" TEMP_DATA_FILENAME = "temp_data.npy" TEMP_SURF_FILENAME = "temp_surf.npy" if len(sys.argv) < 4: sys.exit("Usage: objlake.py npyobjfilename objno depth") print("npyobjfilename=", sys.argv[1]) if os.path.isfile(NPY_DIR+"/"+sys.argv[1]) == False: sys.exit("Object File not found.") m = re.search('npyobj\(([0-9]+)\)\-([0-9]+)\-([0-9]+)\-([0-9]+)\.npy', sys.argv[1]) if m is None: sys.exit("Bad filename") print(m.group(1)) print(m.group(2)) print(m.group(3)) print(m.group(4)) num_obj = int(m.group(1)) mesh1 = int(m.group(2)) mesh2 = int(m.group(3)) mesh3 = int(m.group(4)) objno = int(sys.argv[2]) if objno > num_obj: sys.exit("Bad objno") depth = int(sys.argv[3]) if depth == 0: sys.exit("Bad depth") if os.path.isfile(NPY_DIR+"/"+DATA_FILENAME.format(mesh1, mesh2, mesh3)) == False: sys.exit("Data File not found.") if os.path.isfile(NPY_DIR+"/"+SURF_FILENAME.format(mesh1, mesh2, mesh3)) == False: sys.exit("Surface File not found.") print("Go!") obj_array = np.load(NPY_DIR+"/"+sys.argv[1]) data_array = np.load(NPY_DIR+"/"+DATA_FILENAME.format(mesh1, mesh2, mesh3)) surf_array = np.load(NPY_DIR+"/"+SURF_FILENAME.format(mesh1, mesh2, mesh3)) border_val = data_array.max() + 10 obj_cell = 0 ymin = obj_array.shape[0] ymax = 0 xmin = obj_array.shape[1] xmax = 0 pos1st_y = -1 pos1st_x = -1 print("obj_cell=", obj_cell) print("ymin=", ymin, ", ymax=", ymax) print("xmin=", xmin, ", xmax=", xmax) print("pos1st_y=", pos1st_y, "pos1st_x=", pos1st_x) for y in range(obj_array.shape[0]): for x in range(obj_array.shape[1]): if obj_array[y][x] == objno: obj_cell += 1 if pos1st_y == -1: pos1st_y = y pos1st_x = x if y > ymax: ymax = y if y < ymin: ymin = y if x > xmax: xmax = x if x < xmin: xmin = x print("obj_cell=", obj_cell) print("ymin=", ymin, ", ymax=", ymax) print("xmin=", xmin, ", xmax=", xmax) print("pos1st_y=", pos1st_y, "pos1st_x=", pos1st_x) surf_val = surf_array.max() + 2 data_min = data_array.max() for y in range(ymin, ymax+1): for x in range(xmin, xmax+1): if obj_array[y][x] == objno: if surf_array[y-1][x] == 5: # 上 5=地表面 if data_min > data_array[y-1][x]: data_min = data_array[y-1][x] surf_array[y-1][x] = surf_val if surf_array[y+1][x] == 5: # 下 5=地表面 if data_min > data_array[y+1][x]: data_min = data_array[y+1][x] surf_array[y+1][x] = surf_val if surf_array[y][x-1] == 5: # 左 5=地表面 if data_min > data_array[y][x-1]: data_min = data_array[y][x-1] surf_array[y][x-1] = surf_val if surf_array[y][x+1] == 5: # 左 5=地表面 if data_min > data_array[y][x+1]: data_min = data_array[y][x+1] surf_array[y][x+1] = surf_val print("data_min=", data_min) bottom_val = data_min - depth print("bottom_val=", bottom_val) for y in range(ymin, ymax+1): for x in range(xmin, xmax+1): if obj_array[y][x] == objno: data_array[y][x] = bottom_val np.save(NPY_DIR+"/"+TEMP_LAKE_FILENAME, data_array) if obj_cell > 0: data_array[ymin-2:ymax+2, xmin-2:xmin-1] = border_val data_array[ymin-2:ymax+2, xmax+2:xmax+3] = border_val data_array[ymin-2:ymin-1, xmin-2:xmax+3] = border_val data_array[ymax+2:ymax+3, xmin-2:xmax+3] = border_val np.save(NPY_DIR+"/"+TEMP_DATA_FILENAME, data_array) np.save(NPY_DIR+"/"+TEMP_SURF_FILENAME, surf_array)
後は、出力されたtemp_lake.npyファイルを5x5倍して、池の部分(100:300, 0:300)を切り出してマイクラでロードして確認する。
デバッグ情報として、temp_surf.npyは池の縁の標高として使用した領域を示すデータとして、temp_data.npyは処理した池の範囲を示すboundを2ピクセル大きく表しているデータとして出力しています。
