Skip to content

Commit 29d031b

Browse files
ssmahtossmahto
committed
new commit
1 parent 70f6cb5 commit 29d031b

File tree

7 files changed

+413
-108
lines changed

7 files changed

+413
-108
lines changed

‎code/CURVE.py‎

Lines changed: 51 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -38,61 +38,79 @@ def pick(c, r, mask): # (c_number, r_number, an array of 1 amd 0)
3838
fill.add(south)
3939
return picked
4040

41-
def curve(res_name, max_wl, point, boundary, dem_file_path):
41+
def res_iso(res_name, max_wl, point, boundary, res_directory):
4242
# ===================================================== INPUT PARAMETERS
43-
bc = boundary
44-
pc = point
45-
coords = bc + pc
46-
utm_coords = np.array([utm.from_latlon(coords[i + 1], coords[i]) for i in range(0, len(coords), 2)])
47-
# Bounding box of the reservoir [ulx, uly, lrx, lry]
48-
bbox = np.array([utm_coords[0,0], utm_coords[0,1], utm_coords[1,0], utm_coords[1,1]], dtype=np.float32)
49-
res_point = np.array([utm_coords[2,0], utm_coords[2,1]], dtype=np.float32)
50-
# 30m equivalent distance in degree
51-
#dist30m= 0.01745329252
52-
xp = round(abs(res_point[0]-bbox[0])/30)
53-
yp = round(abs(res_point[1]-bbox[1])/30)
54-
# Maximum reservoir water level
55-
curve_ext = max_wl + 10 # to expand the curve
43+
os.chdir(res_directory + "/Outputs")
44+
res_dem_file = (res_name + "DEM_UTM_CLIP.tif")
45+
dem_ds = gdal.Open(res_dem_file)
46+
geotransform = dem_ds.GetGeoTransform()
5647

57-
# CREATING E-A-S RELATIONSHOP
58-
# clipping DEM by the bounding box
59-
print("Clipping DEM by the bounding box ...")
60-
dem = gdal.Open(dem_file_path)
48+
# Calculate the bounding box coordinates
49+
left = geotransform[0]
50+
top = geotransform[3]
51+
right = left + geotransform[1] * dem_ds.RasterXSize
52+
bottom = top + geotransform[5] * dem_ds.RasterYSize
6153

62-
# Changing path to the desired reservoir
63-
os.chdir(os.getcwd() + "/Outputs")
64-
res_dem_file = res_name+"DEM.tif"
65-
dem = gdal.Translate(res_dem_file, dem, projWin = bbox) #bbox=bc
66-
dem = None
54+
bbox = [left, top, right, bottom]
6755

56+
# 30m nearly equal to 0.00027777778 decimal degree
57+
xp = abs(round((point[0]-boundary[0])/0.00027777778))
58+
yp = abs(round((point[1]-boundary[1])/0.00027777778))
59+
dem_ds = None
60+
61+
# CREATING E-A-S RELATIONSHOP
6862
# isolating the reservoir
6963
dem_bin = gdal_array.LoadFile(res_dem_file)
64+
dem_bin[dem_bin == 32767] = np.nan
7065
#------------------ Visualization <Start>
7166
plt.figure()
7267
plt.imshow(dem_bin, cmap='jet')
7368
plt.colorbar()
74-
dem_bin[np.where(dem_bin > curve_ext)] = 0
69+
70+
dem_bin[np.where(dem_bin > max_wl+10)] = 0 #to expand the reservoir extent for accounting uncertainity in max_wl
7571
dem_bin[np.where(dem_bin > 0)] = 1
7672
res_iso = pick(xp, yp, dem_bin)
73+
plt.figure()
74+
plt.imshow(res_iso, cmap='jet')
75+
plt.colorbar()
7776
output = gdal_array.SaveArray(res_iso.astype(gdal_array.numpy.float32),
7877
"res_iso.tif", format="GTiff",
7978
prototype = res_dem_file)
8079
output = None
8180

8281
# finding the lowest DEM value in the reservoir extent
8382
res_dem = gdal_array.LoadFile(res_dem_file)
84-
res_dem[np.where(res_iso == 0)] = 9999
85-
min_dem = np.nanmin(res_dem)
83+
res_dem[res_dem == 32767] = np.nan
84+
res_dem[np.where(res_iso == 0)] = 9999 # 9999 is any arbitrary unrealistice value
85+
min_dem = int(np.nanmin(res_dem))
8686

87-
# caculating reservoir surface area and storage volume
88-
# coresponding to each water level
87+
88+
#============================================================== E-A relationship
89+
def curve(res_name, res_directory):
90+
# caculating reservoir surface area and storage volume coresponding to each water level
91+
os.chdir(res_directory + "/Outputs")
92+
res_dem_file = (res_name + "DEM_UTM_CLIP.tif")
93+
res_dem = gdal_array.LoadFile(res_dem_file)
94+
res_dem[res_dem == 32767] = np.nan
95+
exp_mask = gdal_array.LoadFile("Expanded_Mask.tif").astype(np.float32)
96+
res_dem[np.where(exp_mask == 0)] = np.nan
97+
output = gdal_array.SaveArray(res_dem.astype(gdal_array.numpy.float32),
98+
"DEM_Landsat_res_iso.tif",
99+
format="GTiff", prototype = res_dem_file)
100+
output = None
101+
# plt.figure()
102+
# plt.imshow(res_dem, cmap='jet')
103+
# plt.colorbar()
104+
min_dem = int(np.nanmin(res_dem))
105+
curve_ext = int(np.nanmax(res_dem)) + 10 # to expand the curve
106+
res_dem_updated = ("DEM_Landsat_res_iso.tif")
107+
89108
results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]
90109
pre_area = 0
91110
tot_stor = 0
92111
for i in range(min_dem, curve_ext):
93112
level = i
94-
water_px = gdal_array.LoadFile(res_dem_file)
95-
water_px[np.where(res_iso == 0)] = 9999
113+
water_px = gdal_array.LoadFile(res_dem_updated)
96114
water_px[np.where(res_dem > i)] = 0
97115
water_px[np.where(water_px > 0)] = 1
98116
area = np.nansum(water_px)*9/10000
@@ -107,8 +125,7 @@ def curve(res_name, max_wl, point, boundary, dem_file_path):
107125
csvWriter = csv.writer(my_csv)
108126
csvWriter.writerows(results)
109127

110-
# ==================== Plot the DEM-based Level-Storage curve
111-
128+
# ==================== Plot the DEM-based Level-Storage curve
112129
data = results[1:, :]
113130
data = np.array(data, dtype=np.float32)
114131
# Extract column 2 and 3 from the array
@@ -124,6 +141,8 @@ def curve(res_name, max_wl, point, boundary, dem_file_path):
124141
plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
125142

126143
return min_dem
144+
145+
127146

128147

129148

‎code/InfeRes.py‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
############+++++++++++ PLEASE MAKE SURE OF THE FOLLOWING POINTS BEFORE RUNNING THE CODE +++++++++++############
2+
# [1]. Data are already downloaded (Satellite images and DEM)
3+
# [2]. DEM should be in the projected coordinate system (unit: meters)
4+
# [3]. Use the same coordinates that you have used in "data_download.py"
5+
# [4]. All the python(scripts) files are inside ".../ReservoirExtraction/codes"
6+
# [5]. "Number_of_tiles" = Number of Landsat tiles to cover the entire reservoir. It is recommended to download the Landsat tile covering the maximum reservoir area
7+
############++++++++++++++++++++++++++++++++++++++++ END ++++++++++++++++++++++++++++++++++++++++#############
8+
9+
# IMPORTING LIBRARY
10+
11+
import os
12+
os.chdir("H:/My Drive/NUSproject/ReservoirExtraction/")
13+
from codes.CURVE import curve
14+
from codes.CURVE import res_iso
15+
from codes.MASK import mask
16+
from codes.WSA import wsa
17+
from codes.PREPROCESING import preprocessing
18+
from codes.CURVE_Tile import one_tile
19+
from codes.CURVE_Tile import two_tile
20+
21+
22+
if __name__ == "__main__":
23+
24+
#====================================>> USER INPUT PARAMETERS
25+
parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
26+
os.chdir(parent_directory)
27+
res_name = "PleiKrong" # Name of the reservoir
28+
res_directory = parent_directory + res_name
29+
# A point within the reservoir [longitude, latitude]
30+
point = [107.872, 14.422]
31+
# Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
32+
boundary = [107.763, 14.672, 107.924, 14.392] #[107.763, 14.672, 107.924, 14.392]
33+
max_wl = 575
34+
yearOFcommission = 2008
35+
Number_of_tiles = 1
36+
os.makedirs(res_name, exist_ok=True)
37+
os.chdir(parent_directory + res_name)
38+
# Create a new folder within the working directory to download the data
39+
os.makedirs("Outputs", exist_ok=True)
40+
# Path to DEM (SouthEastAsia_DEM30m.tif), PCS: WGS1984
41+
dem_file_path = "H:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif"
42+
43+
#====================================>> FUNCTION CALLING -1
44+
# [1]. Data pre-processing (reprojection and clipping)
45+
preprocessing(res_name, point, boundary, parent_directory, dem_file_path)
46+
47+
#====================================>> FUNCTION CALLING -2
48+
# [2]. DEM-based reservoir isolation
49+
res_iso(res_name, max_wl, point, boundary, res_directory)
50+
51+
#====================================>> FUNCTION CALLING -3
52+
# [3]. Creating mask/intermediate files
53+
mask(res_name, yearOFcommission, max_wl, point, boundary, res_directory)
54+
55+
#====================================>> FUNCTION CALLING -4
56+
# [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
57+
res_minElev = curve(res_name, res_directory)
58+
59+
#====================================>> FUNCTION CALLING -5
60+
# [5]. Calculating the water surface area
61+
res_directory = parent_directory + res_name
62+
os.chdir(res_directory)
63+
wsa(res_name, res_directory)
64+
65+
#====================================>> FUNCTION CALLING -6
66+
# [6]. Calculating the reservoir restorage (1 tiles)
67+
if Number_of_tiles==1:
68+
res_directory = parent_directory + res_name
69+
os.chdir(res_directory)
70+
print("One tile reservoir")
71+
one_tile(res_name, max_wl, res_minElev, res_directory)
72+
73+
# Calculation of water surface area for the complete reservoir (2 tiles) and corresponding reservoir restorage
74+
if Number_of_tiles==2:
75+
res_directory = parent_directory + res_name
76+
os.chdir(res_directory)
77+
print("Two tiles reservoir")
78+
# Upper-Left and Lower-right coordinates of the complete reservoir
79+
complete_res_boundary = [101.693, 20.007, 102.331, 19.240]
80+
two_tile(res_name, max_wl, res_minElev, point, complete_res_boundary, dem_file_path, res_directory)
81+
82+
83+
84+
85+
86+
87+
88+
89+
90+
91+
92+
93+
94+

0 commit comments

Comments
 (0)