Python --- Virtually Creating Dams on DEM opend by GSI Japan -Remind memo
Before I will explain the codes and consdertion, I show one of the results that are analysised my codes. The below image is the result that I run the virtual dam creation code near Yamatosaka Dam, 山鳥坂ダム, construction site which is under construction by MILT, Japan at Kawabegawa Branch River Hiji River Ehime Pref.
By Shikoku Regional Development Bureau, MLIT, Yamatosaka Dam's specification is below.(URL: https://www.skr.mlit.go.jp/yamatosa/pamphlet/pdf/yamadam_jigyou.pdf)
- Dam Height: Approximately 96m, 82m (I think that this difference is caused including surface soft soil thickness on the dam site under calculation.)
- Dam Top Length: Approximately 275m, 206m
- Dam Top Elevation: EL161m, EL161m (I try to mutch this values)
- Storage Area: 0.7km2, 669754m2=0.67km2
- Total Storage Volume: 22million cubic meter, 20million cubic meter
It's allowable range for me as rough calculation.
I seem that these conseptis useful to to guess the seservoire volumes which will be made embankment( or natural dams) on river caused by land slide's and the others natural disaster like it. I think that Apps like this have existed in the world.
That aside, I introduce my codes besed on this concept.
At first, we need to import some libraries which are needed our codes.
!pip install geopandas
import folium
from folium import plugins
import sys
from folium.plugins import MarkerCluster
import geopandas as gpd
from shapely.geometry import Point, LineString, MultiPoint, MultiLineString, Polygon
import pandas as pd
At second step, I write the code which have previously introduced the code for downloading Digital Elevation Model of an area from Geospetial Information Authority in Japan.
We need to describe the download map area with north-west corner coordinate as longitude and latitude and necessary number of tiles for direction to south( coodination y) and east(coordination x) like below.
###################################
##### 使用するタイルデータ入力箇所 #####
##### Code No.1 use this data #####
###################################
# Yamatosaka Dam Point
zoom = 13 #zoomを入力
north_west_x = 7114 #北西端tileのx座標を入力
north_west_y = 3283 #北西端tileのy座標を入力
tile_num_x = 6 #東西方向のtile数を入力
tile_num_y = 7 #南北方向のtile数を入力
I pick up Above range of area from GSI Japan's map. It's very wide range of area.
I explain before now about zoom and tile there(URL: https://shikuuk.blogspot.com/2024/10/catchment-area-calculator-with-python.html)(I'm sorry If I mistake the link setting)
I previously made and showed the code to download DEM tiles and to merge tiles
### 国土地理院HPから必要なDEM(標高データ)を取得
### 取得した標高データを図化
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from urllib.request import urlopen, URLError,HTTPError
# 個別タイルのDEM(標高)データを取得する関数'fetch_tile()'
# z:zoom, x:タイルのx座標, y:タイルのy座標
def fetch_tile(z, x, y):
# 変数urlに国土地理院のHPからDEM(標高)データの場所を入力
url = "https://cyberjapandata.gsi.go.jp/xyz/dem/{z}/{x}/{y}.txt".format(z=z, x=x, y=y)
# 指定したURL,FileNameが存在する場合は、
# 変数dfにpandasDataFrameとしてDEMデータを入力し、かつ、
# 標高データに'e'が存在する場合は'0'に置き換え
try:
df = pd.read_csv(url, header=None).replace('e', 0.0)
# 指定したURLやFileNameが存在しない場合（例えば、指定したタイル座標が全て海上の場合）
# 変数dfに全て'0.0'を入力(DEMデータは東西256×南北256のメッシュで構成されており、それ
# らに全て'0.0'を入力)して、海上のタイルに仮のDEMデータを作成
except URLError as e:
temp_tile = []
temp_row = []
zero = 0.0
for i in range(0, 256):
temp_row.append(zero)
for i in range(0, 256):
temp_tile.append(temp_row)
df = pd.DataFrame(temp_tile)
except HTTPError as e:
temp_tile = []
temp_row = []
zero = 0.0
for i in range(0, 256):
temp_row.append(zero)
for i in range(0, 256):
temp_tile.append(temp_row)
df = pd.DataFrame(temp_tile)
return df.values
# 北西端・南東端タイル座標を指定した矩形領域の複数タイルのDEM(標高)データを
# 結合する関数'fetch_all_tiles'
def fetch_all_tiles(north_west, south_east):
assert north_west[0] == south_east[0] , "タイル座標のzが一致していません"
x_range = range(north_west[1], south_east[1] + 1)
y_range = range(north_west[2], south_east[2] + 1)
return np.concatenate(
[
np.concatenate(
[fetch_tile(north_west[0], x, y) for y in y_range],
axis = 0
) for x in x_range
],
axis = 1
)
# 南東端のタイル番号の計算
south_east_x = north_west_x + tile_num_x -1 #南東端tileのx座標
south_east_y = north_west_y + tile_num_y -1 #南東端tileのy座標
# 複数タイルの統合
tile = fetch_all_tiles((zoom, north_west_x, north_west_y), (zoom, south_east_x, south_east_y))
tile = tile.astype(float) #str型の数値をfloat型に変換
We try to show merged DEM tiles to the plot of the matplotlib graph by below code.
### 取得した標高データを図化
# 等高線の作図
# 等高線のx,y軸グリッドを作成
grid_x = tile_num_x * 256
grid_y = tile_num_y * 256
X, Y = np.meshgrid(np.linspace(grid_x, 0, grid_x), np.linspace(grid_y, 0, grid_y))
fig = plt.figure(figsize=(14,10), dpi=80, facecolor = 'w', edgecolor = 'k')
plt.gca().invert_yaxis()
# 等高線を描く
elevation = range(0,2200,50) #contourのmin/maxとpitchを指定
cont = plt.contour(tile, levels = elevation, cmap='cool')
#cont = plt.contour(tile, cmap='spring_r') # c_mapの種類の違い
# 等高線の凡例をつける
cb = plt.colorbar(cont, shrink=0.5, aspect=10)
plt.show()
Next step, I define some function to operate coordination(examply transformation from pixel coordinates to longitude, latitude and reversive exchange and the others. The transformation between pixel coordinate of GSI Japan's tile and Longitude, Latitude has complex equation, previously told.( ex. https://shikuuk.blogspot.com/2025/02/pythondemcodefunction.html)
### code No2 ###
### 経度(longitude)・緯度(latitude)とLocalTileCoordinate間の座標変換
### タイル座標からタイル北西端の緯度経度を計算
### タイル内部座標のメッシュサイズなどの計算
### タイル標高値の平滑化処理
!pip install pyproj
from math import pi
from math import e
from math import atan
import pyproj
# Tile北西端のGlobalTileCoordinate(PixelCoord)から緯度Lat経度LonをCalculateするFunction
# z:zoom, x:GlobalTileCoord(PixelCoord)のx座標, y:GlobalTileCoord(PixelCoord)のy座標
def func_tile_to_latlon(z, x, y):
lon = (x / 2.0**z) * 360 - 180 # Lon（東経）
mapy = (y / 2.0**z) * 2 * pi - pi
lat = 2 * atan(e ** (- mapy)) * 180 / pi - 90 # Lat（北緯）
return lon,lat
# LocalTile内のMeshSize、MeshArea(Unit:m2)をCalculateするFunction
##### このCodeのRequirementはのちほど精査 #####
def func_mesh_size_inner_tile(z, x, y):
lon1, lat1 = func_tile_to_latlon(z, x, y)
lon2, lat2 = func_tile_to_latlon(z, x+1, y+1)
lon_size = (lon2-lon1)/256
lat_size = (lat2-lat1)/256
lon_pich = func_lonlat_to_dist(lon1, lat1, lon2, lat1) /256
lat_pich = func_lonlat_to_dist(lon1, lat1, lon1, lat2) /256
area = lon_pich * lat_pich
return lon_size, lat_size, lon_pich, lat_pich, area
# 任意のLocalColumnRowNumber'cr'からLonLat(経度緯度)をCalculateするFunction'func_colrow_to_lonlat()'
def func_colrow_to_lonlat(zoom, colrow, nw_edge_lonlat):
parameter_l = 85.05112878
col, row = colrow
nw_edge_lon, nw_edge_lat = nw_edge_lonlat
# LocalTile北西端のPixelCoordをCalculate
nw_edge_col = int(2 ** (zoom + 7) * (nw_edge_lon / 180 + 1))
nw_edge_row = int(2 ** (zoom + 7) / pi * (-1 * np.arctanh(np.sin(pi * nw_edge_lat / 180)) + np.arctanh(np.sin(pi / 180 * parameter_l))))
# LocalTileのColumnRowNumberをPixelCoordにExchange
pixcel_col = col + nw_edge_col
pixcel_row = row + nw_edge_row
lon = 180 * (pixcel_col / 2 ** (zoom + 7) - 1)
lat = 180/pi * (np.arcsin(np.tanh(-1 * pi / (2**(zoom+7)) * pixcel_row + np.arctanh(np.sin(pi/180 *parameter_l)))))
return lon, lat
# Function that exchange to column and row number in the conbined elevation model from zoom and north-west edge lon,kat
#/ 任意の緯度・経度について、結合した標高値モデルのzoomと北西端緯度・経度から結合標高値モデルにおける行列番号に変換する関数
#/ Latest version as of 20251121
def func_lonlat_to_colrow(z, lonlat, nw_edge_lonlat):
parameter_l = 85.05112878
lon, lat = lonlat
wn_edge_lon, wn_edge_lat = nw_edge_lonlat
# Calculate inner(local) tile column, row and the Pixel Coordinate(Global Coordinate) from zoom and lon, lat
# / zoomと緯度経度から全地球系のPixcel座標を計算
global_column = int(2 ** (z + 7) * (lon / 180 + 1))
global_row = int(2 ** (z + 7) / pi * (-1 * np.arctanh(np.sin(pi * lat / 180))+np.arctanh(np.sin(pi / 180 * parameter_l))))
# Calculate the Pixel coordinate of west and north edge of the combined elevation model
w_edge_column = int(2 ** (z + 7) * (wn_edge_lon / 180 + 1))
n_edge_row = int(2 ** (z + 7) / pi * (-1 * np.arctanh(np.sin(pi * wn_edge_lat / 180)) + np.arctanh(np.sin(pi / 180 * parameter_l))))
# Calculate the row and column number in the local coodinate system(at conbined elevation model)
column = global_column - w_edge_column
row = global_row - n_edge_row
return column ,row
# LotLan(経度緯度)が与えられた2点間DistanceをCalculateするFunction'func_lonlat_to_dist()'
def func_lonlat_to_dist(lon1, lat1, lon2, lat2):
grs80 = pyproj.Geod(ellps='GRS80')
az, bkw_az, dist = grs80.inv(lon1, lat1, lon2, lat2)
return dist
####### 対象領域の北西・南東端の経度・緯度の計算
lonlat_nw_edge = func_tile_to_latlon(zoom, north_west_x, north_west_y)
lonlat_se_edge = func_tile_to_latlon(zoom, south_east_x + 1, south_east_y + 1)
###### タイル内部のメッシュサイズの計算 ########
lon_mesh_size, lat_mesh_size, lon_mesh_pich, lat_mesh_pich, mesh_area = func_mesh_size_inner_tile(zoom, north_west_x, north_west_y)
# 計算確認用出力
# タイルに関する諸元の出力
print("タイル北西端経度緯度", lonlat_nw_edge, "北西端タイル番号x y",north_west_x, north_west_y)
print("mesh_size", lon_mesh_size, lat_mesh_size)
print("mesh_area:m2", mesh_area)
print('Some essential function has been run')
the output of above code is below.
Requirement already satisfied: pyproj in /usr/local/lib/python3.12/dist-packages (3.7.2)
Requirement already satisfied: certifi in /usr/local/lib/python3.12/dist-packages (from pyproj) (2025.11.12)
タイル北西端経度緯度 (132.626953125, 33.61461929233377) 北西端タイル番号x y 7114 3283
mesh_size 0.000171661376953125 -0.00014298650855371076
mesh_area:m2 252.64218331784716
Some essential function has been run
Until here, We can take some basic codes for operation the GeoSpacial Information Authority of Japan.
I'm finally starting vertuial dam creation on the GSI DEM map.
At first on this step, we should make a code to create some linier lines with a given point as the center and changing the slope by 0.2 and set some wall for upper and lower side each line as the initial point to calculate the reservoir area.
Meanwhile Extend each line in a certain direction, and set the point where the elevation of the extension exceeds the dam crest elevation, which is the elevation of the point plus the dam embankment height, as the short point of the line. And distances each lines' two short points become the vurtual dam top lengthes.
def func_temp_dam_walls_create(col_0, row_0):
global dam_wall_cr_list #crはcolrowの略語
global ini_sa_cr_list1
global ini_sa_cr_list2
global temp_ini_sa_cr_line_list
global bound_line_cr_lists
###########################################
# Tentitive Dam Wall & #
# Initial Strage Calculate Points Created #
###########################################
# DamWallのLocalTileCol,Rowを収納するdam_wall_cr_listを初期化
dam_wall_cr_list = []; ini_sa_cr_list1 = []; ini_sa_cr_list2 = [] ; temp_ini_sa_cr_line_list = []
for i in range(-5,5):
line_coord_cr_list = []
ini_sa_cr_list1 = []; ini_sa_cr_list2 = []
for m in range(-1, 2, 2):
col = col_0; row = row_0
length = 0.0
dam_lonlat = func_colrow_to_lonlat(zoom, [col_0, row_0], lonlat_nw_edge)
while length <= limit_dam_length and tile[int(row), int(col)] <= dam_top_elev:
""" ダム地点(col_0, row_0)と現在確認中のpixel座標(col, row)の距離lengthが
初期設定した最大長limit_dam_lengthより小さい
または、
現在確認中のpixel座標(col, row)の標高tile[int(row), int(col)]が
初期設定したダム堤頂標高dam_top_elevより小さければ以下のCodeをRun """
colrow = [int(col), int(row)]
lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
length = func_lonlat_to_dist(lon, lat, dam_lonlat[0], dam_lonlat[1])
for j in range(-2, 3):
""" j==-2、j==2の場合は、Dam_Wallに隣接する上下流側に貯水池評価
の起点となる座標としてini_sa_cr_list1,ini_sa_cr_list2に格納
この時点でlist1 or list2のどちらがDamの上流側か、下流側かは未定
J=-1, 0, 1の場合は、Dam本体(Dam_wall)を構成するpixel座標として
line_coord_cr_listにPixel座標を格納 """
# Storage_Areaの評価起点ini_sa_cr_list1とini_sa_cr_list2を設定
if j == -2:
ini_sa_cr_list1.append([col, row + j])
elif j == 2:
ini_sa_cr_list2.append([col, row + j])
# Dam_Wallを設定
else:
line_coord_cr_list.append([col, row + j])
""" 初期座標(col_0, row_0)を中心を通る勾配i*0.2の直線に沿って、現在の
Pixel座標に隣接する座標点を計算 """
col = col + m
row = int(row_0 + (col - col_0) * (i * 0.2))
print(f"row= {int(row)}, col= {int(col)}")
""" 初期座標(col_0, row_0)、列に対して勾配i*0.2のDam_wallを構成する座標
を格納"""
dam_wall_cr_list.append(line_coord_cr_list)
temp_ini_sa_cr_line_list.append([ini_sa_cr_list1, ini_sa_cr_list2])
""" 先ほどの同様のOperationを行に対して勾配i*0.2で計算 """
for i in range(-5,5):
line_coord_cr_list = []
ini_sa_cr_list1 = []; ini_sa_cr_list2 = []
for m in range(-1, 2, 2):
col = col_0; row = row_0
length = 0.0
while length <= limit_dam_length and tile[int(row), int(col)] <= dam_top_elev:
colrow = [int(col), int(row)]
lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
length = func_lonlat_to_dist(lon, lat, dam_lonlat[0], dam_lonlat[1])
for j in range(-2, 3):
# Dam_Wallの上下流に貯水池評価の起点となる座標を設定
if j == -2:
ini_sa_cr_list1.append([col + j, row])
elif j == 2:
ini_sa_cr_list2.append([col + j, row])
# Dam_Wallを設定
else:
line_coord_cr_list.append([col + j, row])
row = row + m
col = int(col_0 + (row - row_0) * (i * 0.2))
dam_wall_cr_list.append(line_coord_cr_list) #
temp_ini_sa_cr_line_list.append([ini_sa_cr_list1, ini_sa_cr_list2])
# print(temp_ini_sa_cr_line_list)
# print("length of dam_wall_cr_list = ", len(dam_wall_cr_list))
########################
# BoundaryWall Created #
########################
bound_line_cr_lists = []
for i in range(-5, 5):
line_coord_cr_list = []
for m in range(-1, 2, 2):
col = col_0; row = row_0
length = 0.0
while length <= limit_bound_length and tile[int(row), int(col)] <= dam_top_elev:# + afford_hight:
colrow = [int(col), int(row)]
lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
length = func_lonlat_to_dist(lon, lat, dam_lonlat[0], dam_lonlat[1])
for j in range(-1, 2):
## Dam_Wallの上下流に貯水池評価の起点となる座標を設定
line_coord_cr_list.append([col, row + j])
col = col + m
row = int(row_0 + (col - col_0) * (i * 0.2))
bound_line_cr_lists.append(line_coord_cr_list)
for i in range(-5, 5):
line_coord_cr_list = []
ini_sa_cr_lists = []
for m in range(-1, 2, 2):
col = col_0; row = row_0
length = 0.0
while length <= limit_bound_length and tile[int(row), int(col)] <= dam_top_elev:# + afford_hight:
colrow = [int(col), int(row)]
lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
length = func_lonlat_to_dist(lon, lat, dam_lonlat[0], dam_lonlat[1])
for j in range(-1, 2):
# Dam_Wallの上下流に貯水池評価の起点となる座標を設定
# Dam_Wallを設定
line_coord_cr_list.append([col + j, row])
row = row + m
col = int(col_0 + (row - row_0) * (i * 0.2))
bound_line_cr_lists.append(line_coord_cr_list) #
return
The points of innovation in this code is to set the 'ini_sa_cr_list1' and 'ini_sa_cr_list2' at the upper and lower side of the vartual dam wall, to use either list of coordinates as the starting point for calculating the resevoir's(straging) planer extent.
I'm afraid to say that I forget what difference between 'temp_ini_sa_cr_line_list' and 'bound_line_cr_lists'.(Maybe the boundary is a line that should not be crossed when calculating the reservoir area.)
Next step we should select the dam that has minimum dam top length, from dam linears has a various of direction on the 2D plane which previously define as dam walls.
This means that Select the shortest line segment from among multiple virtual dam line segments drawn radially from a certain point.
'min_dam_length = limit_dam_length * 1.5' In this code, means that we give min_dam_length limit_dam_length * 1.5, where limit_dam_length is the threshold value setting vertual dam_walls. if dam_top_length exceed this value, the relevant dam_wall is omited from the consider, so that dam_top_length too long.
'dam_wall_cr_list' is the pixel coordinate values list of multiple dam wall line segments created radially from a given point are stored.
I use the 'for' loop to select the dam wall which have shortest length of the dam_top_length. At this time, I use 'itertools.combinations()' function that it's very powerful and convenience tool.
If we want to know the details, we try to ask Master Gemini by copy and paste the below codes.
#Dam Top LengthがMinimumなDam Wall(Line)をSelectするFunction
import itertools
def select_min_len_Dam_Wall():
global dam_body_cr_list
global ini_sa_cr_lists
global bound_line_cr_list
global min_dam_length
#最小堤長を仮設定
min_dam_length = limit_dam_length * 1.5
#最小堤長となるColRowListがdam_wall_cr_listの何番目にあたるかを設定するための仮整数
n = 0; num_of_min_line = 0
#堤長が最も短いDamLineを設定
for dam_line_cr_list in dam_wall_cr_list:
#最初の2点を仮の最大距離点として設定
colrow1, colrow2 = dam_line_cr_list[0], dam_line_cr_list[1]
lon1, lat1 = func_colrow_to_lonlat(zoom, colrow1, lonlat_nw_edge)
lon2, lat2 = func_colrow_to_lonlat(zoom, colrow2, lonlat_nw_edge)
line_length = func_lonlat_to_dist(lon1, lat1, lon2, lat2)
#すべての座標の組み合わせに対して距離を計算
for colrows in itertools.combinations(dam_line_cr_list, 2):
lon1, lat1 = func_colrow_to_lonlat(zoom, colrows[0], lonlat_nw_edge)
lon2, lat2 = func_colrow_to_lonlat(zoom, colrows[1], lonlat_nw_edge)
length = func_lonlat_to_dist(lon1, lat1, lon2, lat2)
if length > line_length:
line_length = length
# print("line_length", line_length)
if line_length < min_dam_length:
min_dam_length = line_length
dam_body_cr_list = dam_line_cr_list.copy()
# print("ダム本体の堤長を", line_length, "に置き換えました")
num_of_min_line = n
n = n + 1
# print("ダム本体の堤長が最小となるのは")
# print(dam_body_cr_list)
# print("ダム堤長は", min_dam_length, "mです。")
# print("length of dam_body_cr_list = ", len(dam_body_cr_list))
ini_sa_cr_lists = temp_ini_sa_cr_line_list[num_of_min_line]
bound_line_cr_list = bound_line_cr_lists[num_of_min_line]
# dam_body_lonlat_list = []
# for colrow in dam_body_cr_list:
# lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
# dam_body_lonlat_list.append([lon, lat])
# bound_line_lonlat_list = []
# for colrow in bound_line_cr_list:
# lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
# bound_line_lonlat_list.append([lon, lat])
# print("length of ini_sa_cr_list = ", len(ini_sa_cr_lists[1]))
return
Next, I define the function that calculates the dam top length like below. This logic like this have been used above code.
# DamTopLengthをCalするFunction...dam_top_len()
def func_dam_top_len():
global dam_body_cr_list
length = 0; dam_top_length = 0
for colrows in itertools.combinations(dam_body_cr_list, 2):
lon1, lat1 = func_colrow_to_lonlat(zoom, colrows[0], lonlat_nw_edge)
lon2, lat2 = func_colrow_to_lonlat(zoom, colrows[1], lonlat_nw_edge)
length = func_lonlat_to_dist(lon1, lat1, lon2, lat2)
if length > dam_top_length:
dam_top_length = length
return dam_top_length
Next, I try to calculate the resevoir's(water strage) area and volume. Please look at the below three function.
The first function 'func_sa_process(c_num, r_num)’ selects the eight points around a given pixel coordinate(c_num, r_num), and if the elevation of each point is equal to or less than the elevation of the dam, the point is evaluated as a reservoir and the pixel coordinate of the point add to 'sa_cr_eval_list'.
The second fnction 'func_sa_points()' gives the candidate point that has potential for becoming strage area to 'func_sa_process()’ in order, and make tha final storage pixel coordinate list 'sa_cr_list'.
The third function 'func_sa_area_volume_cal(sa_cr_list, dam_top_elev)' calculates the strage(reservoir) area and storage(reservoir) volum from 'sa_cr_list'
If we want to know the detail of these codes and function, we could ask the Master
Gemini.
#任意の内部タイル座標の隣接内部タイル標高が貯水面以下にあることを調べる関数'func_sa_process'
#'sa' = storage_area
def func_sa_process(c_num, r_num):
global sa_cr_eval_list
global temp_sa_cr_list
global sa_cr_list
n = 0
# print(c_num, r_num)
for i in range(-1, 2):
for j in range(-1,2):
n += 1
# 当該タイルは評価しない
# print(tile[r_num + i][c_num + j], dam_top_elev)
if i == 0 and j == 0:
continue
## 比較対象タイルがダウンロードしたタイル範囲を超えたら評価しない
# if (c_num + i < 0) or (c_num + i >= x_num):
# continue
# if (r_num + j < 0) or (r_num + j >= y_num):
# continue
if tile[r_num + i][c_num + j] <= dam_top_elev - 0: #########マイナス値はダム堤頂高さに対する余裕高さ
# print("+++col, row", c_num + j, r_num + i, "Ground Level", tile[r_num + i][c_num + j], "dam_ground_level", dam_ground_level)
# sa_elav_listの整理
if [c_num + j, r_num + i] in sa_cr_eval_list:
pass
####20231002delete elif [c_num + i, r_num + j] in dam_body_cr_list:
elif [c_num + j, r_num + i] in bound_line_cr_list:
pass
# 2023/09/19 Adding from here
elif tile[r_num + i][c_num + j] < dam_ground_level - 50 or tile[r_num + i][c_num + j] <= 0.0 \
or r_num <= 0 or c_num <= 0 or r_num >= 256 * tile_num_y or c_num >= 256* tile_num_x :
###########################################################
print("sa_cr_list", sa_cr_list)
# print("++++n = ", n) ####################
sa_cr_eval_list = []
temp_sa_cr_list = []
# print("***EEEE***", sa_cr_eval_list)
#######20230927_Caution!!!! temp_sa_cr_list = [] #20230922_adding
return
# 2023/09/19 Adding to here
else:
sa_cr_eval_list.append([c_num + j, r_num + i])
# ca_eval_numから隣接内部タイルとの高低差確認を実施したタイル内部座標(c_num, r_num)を削除
if [c_num, r_num] in sa_cr_eval_list:
sa_cr_eval_list.remove([c_num, r_num])
return
# 湛水領域を格納していく関数［global関数のみを使用]
def func_sa_points():
global sa_cr_eval_list
global temp_sa_cr_list
global sa_cr_list
while len(sa_cr_eval_list) != 0:
# print("+++AAA++++", sa_cr_eval_list) ####
for inner_tile_coord in sa_cr_eval_list:
# print("+++PPP+++sa_cr_eval_list", sa_cr_eval_list)
if inner_tile_coord in temp_sa_cr_list:
if inner_tile_coord in sa_cr_eval_list:
sa_cr_eval_list.remove(inner_tile_coord)
continue
elif inner_tile_coord in bound_line_cr_list:
if inner_tile_coord in sa_cr_eval_list:
sa_cr_eval_list.remove(inner_tile_coord)
continue
# 評価する内部タイルの内部座標をキャッチメントリスト'ca_inner_tile_cr_list'に追加
else :
temp_sa_cr_list.append(inner_tile_coord)
func_sa_process(inner_tile_coord[0], inner_tile_coord[1])
if len(sa_cr_eval_list) == 0:
return
# print("+++JJJ+++sa_cr_eval_list", sa_cr_eval_list)
sa_lonlat = func_colrow_to_lonlat(zoom, [inner_tile_coord[0], inner_tile_coord[1]], lonlat_nw_edge)
if len(temp_sa_cr_list) > 0:
sa_cr_list = temp_sa_cr_list.copy()
# print("+FFFFF+", len(sa_cr_eval_list))
# 湛水面積と湛水容量を計算する関数
def func_sa_area_volume_cal(sa_cr_list, dam_top_elev):
area = len(sa_cr_list) * mesh_area
volume = 0.0
for colrow in sa_cr_list:
volume += (dam_top_elev - tile[colrow[1]][colrow[0]]) * mesh_area
# print('####', dam_top_elev - tile[colrow[1]][colrow[0]])
return area, volume # Units are area:m2, volume:m3
But I defined I seem not to be essential function 'transform_colrow_to_lonlat()', but I wrote its function in main code, just in case I introduce the code.
# LocatTileCoordinateをLonLanCoodinateにTransformするFunction
def transform_colrow_to_lonlat():
global sa_lonlat_list
global ini_sa_cr_line_lonlat_list
global dam_body_lonlat_list
global bound_line_lonlat_list
# TestOutput用Code
for coord in sa_cr_list:
sa_lonlat = func_colrow_to_lonlat(zoom, [coord[0], coord[1]], lonlat_nw_edge)
sa_lonlat_list.append([sa_lonlat[0], sa_lonlat[1]])
ini_sa_cr_line_lonlat_list = []
for coord in ini_sa_cr_lists[0]:
sa_lonlat = func_colrow_to_lonlat(zoom, [coord[0], coord[1]], lonlat_nw_edge)
ini_sa_cr_line_lonlat_list.append([sa_lonlat[0], sa_lonlat[1]])
for coord in ini_sa_cr_lists[1]:
sa_lonlat = func_colrow_to_lonlat(zoom, [coord[0], coord[1]], lonlat_nw_edge)
ini_sa_cr_line_lonlat_list.append([sa_lonlat[0], sa_lonlat[1]])
dam_body_lonlat_list = []
for colrow in dam_body_cr_list:
lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
dam_body_lonlat_list.append([lon, lat])
bound_line_lonlat_list = []
for colrow in bound_line_cr_list:
lon, lat = func_colrow_to_lonlat(zoom, colrow, lonlat_nw_edge)
bound_line_lonlat_list.append([lon, lat])
So it's coming near the ending stage.
Next, I write a code of the fuction 'add_result_to_map(dam_lonlat)' that show the result of the dam evaluation and that by follium.
def add_result_to_map(dam_lonlat):
global map
folium.Marker(location=[dam_lonlat[1], dam_lonlat[0]], icon=folium.Icon(color="red"), \
popup = '<i>lon' + str(dam_lonlat[0]) + ', lat ' + str(dam_lonlat[1]) + \
', Dam Height(m)' + str(dam_height) + ', Dam Top Elev.(EL:m)' + str(int(dam_top_elev)) + \
', Dam Top Length(m)' + str(int(dam_top_len)) + ', Strage Area(m2) ' + str(int(sa_area)) + \
', Strage_Volume(m3) ' + str(int(sa_volume)),
).add_to(map)
for data in bound_line_lonlat_list:
folium.Circle(radius=8, location=[data[1], data[0]], color="grey", fill=True, fill_opacity=1).add_to(map)
folium.Circle(radius=8, location=[data[1], data[0]], color="grey", popup = '<i>lon' + str(data[0]) + ', lat ' + str(data[1]), fill=True, fill_opacity=1).add_to(map)
for data in dam_body_lonlat_list:
folium.Circle(radius=8, location=[data[1], data[0]], color="black", fill=True, fill_opacity=1).add_to(map)
for sa in sa_lonlat_list:
folium.Circle(radius=10, location=[sa[1], sa[0]], weight=0, color="blue", fill=True).add_to(map)
for sa in ini_sa_cr_line_lonlat_list:
folium.Circle(radius=10, location=[sa[1], sa[0]], weight=0, color="pink", fill=True).add_to(map)
##地図保存先
#map.save('map.html')
def func_add_eval_result_to_map(initial_lonlat):
""" 計算条件を満たさず、貯水池評価ができなかった場合に評価した初期地点だけ
をFoliumMap上に表示するコード"""
global map
folium.Marker(location=[initial_lonlat[1], initial_lonlat[0]], icon=folium.Icon(color="gray"), \
popup = '<i>lon' + str(initial_lonlat[0]) + ', lat ' + str(initial_lonlat[1]) + \
', Dam Height(m)' + str(dam_height) + ', Dam Top Elev.(EL:m)' + str(int(dam_top_elev)) + \
', Dam Top Length(m)' + str(int(dam_top_len)) + ', Strage Area(m2) ' + str(int(sa_area)) + \
', Strage_Volume(m3) ' + str(int(sa_volume)),
).add_to(map)
I should give the code some initial input data about the considering point logtude, latitude and the folium setting like below.
#######################################
### Input Data for Main Code No.2 ###
#######################################
# Calculate some values for dam point
# dam pointの緯度経度
eval_point_lonlat_list = [[132.70722610796057, 33.4600296142955] #山鳥坂ダム地点
]
# Initialize Folium's Variable'map'
map = folium.Map(location=[eval_point_lonlat_list[0][1], eval_point_lonlat_list[0][0]],
tiles = 'https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png',
attr = '国土地理院 地球地図',
zoom_start = 14)
At final stage, I wrote the main code which connected some functions, I difined before.
########################
#### Main Code No.2 ####
########################
""" Main Code No.2は小さな領域にある
少数地点を評価するためのコードです"""
########################
# Show Initial Setting #
########################
for eval_point_lonlat in eval_point_lonlat_list:
# Calculate col, row some values for point to be evaluated
# Exchange from lon, lat of the point to be evaluated to pixel coordination col, row
col, row = func_lonlat_to_colrow(zoom, eval_point_lonlat, lonlat_nw_edge)
# initial evalatied point to col_0, row_0 as Initial point value
col_0, row_0 = col, row
# Col_0とRow_0のElevationをGet
dam_ground_level = tile[row_0, col_0] # EL:meter
# DamTopLevel(EL:m)をCalculate
dam_top_elev = dam_ground_level + dam_height # EL:meter
print("dam_top_elev", dam_top_elev)
########################
# Main Code below here #
########################
# Initialize a varius of Valuables
dam_wall_cr_list = []; ini_sa_cr_list1 = []; ini_sa_cr_list2 = [] ; temp_ini_sa_cr_line_list = []
bound_line_cr_lists = []; dam_body_cr_list = []; ini_sa_cr_lists = []; bound_line_cr_list = []
sa_cr_list = []; sa_lonlat_list = []; temp_sa_cr_list = []; sa_cr_eval_list = []
print("Longitude, Latitude of the given point", eval_point_lonlat)
print("Pixel Coordinate of the given dam point", int(row), int(col))
func_temp_dam_walls_create(col_0, row_0)
select_min_len_Dam_Wall()
if len(dam_body_cr_list) > 0:
dam_top_len = func_dam_top_len()
else:
func_add_eval_result_to_map(eval_point_lonlat)
pass
if dam_top_len > 2000: # If dam_top_length is more than limit length, the point don't evaluate.
print("Dam top length is over 1000m. Storage Area calculation don't execute!!!")
max_ave_elev = 0; max_ave_elev_ini_sa_cr_list = []
for ini_sa_cr_list in ini_sa_cr_lists:
""" Dam_Wallの両側２つのini_sa_cr_listのうち、平均標高が高い方をSelectして、
貯水池評価計算の初期点となる側のini_sa_cr_listをsa_cr_eval_listにCopy"""
num_of_coord = 0; total_elev = 0
for ini_sa_cr_coord in ini_sa_cr_list:
total_elev += tile[ini_sa_cr_coord[1], ini_sa_cr_coord[0]]
num_of_coord += 1
current_ave_elev = total_elev / num_of_coord
print("current average Elevation(EL:m)", current_ave_elev)
if current_ave_elev > max_ave_elev:
max_ave_elev = current_ave_elev
max_ave_elev_ini_sa_cr_list = ini_sa_cr_list.copy()
sa_cr_eval_list = max_ave_elev_ini_sa_cr_list.copy()
print("++++++++", ini_sa_cr_list)#####
func_sa_points()
transform_colrow_to_lonlat()
sa_area, sa_volume = func_sa_area_volume_cal(sa_cr_list, dam_top_elev)
# dam_body_vol = dam_height / 3 * dam_height / 2 * min_dam_length / 2 # 重力式ダムの場合
""" 重力式ダムの場合のダム堤体積は最大堤高部分の横断面を高さ＝ダム高さ、底辺はダム高さの
３分の１の三角形と仮定し、さらに、ダム端部は断面積をゼロとして、平均断面法で求めた
体積と仮定して計算。"""
dam_body_vol = (dam_height ** 2) / 2 * min_dam_length / 2 # 天然ダムの場合
""" 天然ダムの場合のダム堤体積は最大堤高部分の横断面を高さ＝底辺＝ダム高さの
三角形と仮定し、さらに、ダム端部は断面積をゼロとして、平均断面法で求めた
体積と仮定して計算。"""
# ダムが成立するならば、ダムと貯水池をmapに入力し、
# ダムが成立しないなら、計算した座標のみをmapに入力する。
if len(sa_cr_list) == 0:
func_add_eval_result_to_map(eval_point_lonlat)
else:
add_result_to_map(eval_point_lonlat)
print("Storage(Reservoir) Area Pixel Coordinate List", sa_cr_list)
#地図保存先
map.save("dam_test4.html")
#地図の表示
map
We could get the storage(reservoire)'s map like below.
And I regert to write a code that select which col and row list is the side of reservoire for the given dam wall in the main code. I should write the code as a fuction.
I think that I try to tranlate this code to Web application, so I do the functionilze the code at that time.
At this time, I want to commite the version of this code by Giit.
I run the code below in order to confirm folder's status by !git status.
I have many untracked fill which listed after avove image, but I omit them.
I run to !git add and !git commit -m like below.
But I want to be better code here. Here is the essential judge for calculating the pixel coordinates that incliude the target Storage(Reservoire) Area. But the codes here should not include the main code.
I think of writing this activity to the different BLOG.
URL to Next Page
コメント
コメントを投稿