Skip to content

Commit 9ef4ff6

Browse files
ssmahtossmahto
committed
new_commit
1 parent 16297a2 commit 9ef4ff6

File tree

8 files changed

+257
-110
lines changed

8 files changed

+257
-110
lines changed

‎code/CURVE.py‎

Lines changed: 18 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ def best_fit_degree (xdata, ydata):
5959
from sklearn.metrics import r2_score
6060

6161
# Maximum degree for polynomial regression
62-
max_degree = 5
62+
max_degree = 4
6363

6464
# Lists to store R-squared values and RMSEs for each degree
6565
r_squared_values = []
@@ -154,8 +154,8 @@ def curve_extrapolate(xdata, ydata, bf_degree):
154154
coefficients = np.polyfit(xdata, ydata, degree)
155155
p = np.poly1d(coefficients)
156156

157-
# Generate area values approaching zero for extrapolation
158-
extrapolated_xdata_values = np.linspace(0, max(xdata), 1000)
157+
# Generate area values approaching 30% of the largest area for extrapolation
158+
extrapolated_xdata_values = np.linspace(0.2*min(xdata), max(xdata), 1000)
159159
#extrapolated_area_values = np.linspace(max(area), 1.5*max(area), 1000)
160160
# Use the polynomial function to extrapolate elevation values for small area values
161161
extrapolated_ydata_values = p(extrapolated_xdata_values)
@@ -267,8 +267,8 @@ def curve_preDEM(res_name, point_loc, res_directory):
267267
# Extract column 2 and 3 from the array
268268
data = results[1:, :]
269269
data = np.array(data, dtype=np.float32)
270-
area = data[:, 1] # area
271-
elevation = data[:, 0] # elevation (ydata)
270+
area = data[2:, 1] # area
271+
elevation = data[2:, 0] # elevation (ydata)
272272

273273
# Function calling to get best-fit degree of polinomial
274274
bf_degree = best_fit_degree(area, elevation)
@@ -308,15 +308,15 @@ def curve_preDEM(res_name, point_loc, res_directory):
308308
# Set labels and title
309309
plt.xlabel('Level (m)')
310310
plt.ylabel('Storage (mcm)')
311-
plt.title(res_name)
311+
plt.title(res_name + ' (Minimum DEM level= '+ str(round(data[0,0]))+'m)')
312312
plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
313313

314314
return round(data[0,0])
315315

316316

317317

318318
#============================================================== E-A relationship
319-
def curve_postDEM(res_name, res_directory):
319+
def curve_postDEM(res_name, max_wl, res_directory):
320320
# caculating reservoir surface area and storage volume coresponding to each water level
321321
os.chdir(res_directory + "/Outputs")
322322
res_dem_file = (res_name + "_DEM_UTM_CLIP.tif")
@@ -334,7 +334,7 @@ def curve_postDEM(res_name, res_directory):
334334
# plt.colorbar()
335335
#
336336
min_dem = int(np.nanmin(res_dem))
337-
curve_ext = int(np.nanmax(res_dem))
337+
curve_ext = max_wl+20
338338
res_dem_updated = ("DEM_Landsat_res_iso.tif")
339339

340340
results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]
@@ -355,8 +355,12 @@ def curve_postDEM(res_name, res_directory):
355355
# Extract column 2 and 3 from the array
356356
data = results[1:, :]
357357
data = np.array(data, dtype=np.float32)
358-
area = data[:, 1] # area
359-
elevation = data[:, 0] # elevation (ydata)
358+
idx = np.where(data[:, 1]==min(data[:, 1]))
359+
data = data[np.size(idx,1)-1:]
360+
idx = np.where((data[:, 0] > 0) & (data[:, 1] > 0) & (data[:, 2] > 0))[0]
361+
data = data[idx[0]:]
362+
area = data[2:, 1] # area
363+
elevation = data[2:, 0] # elevation (ydata)
360364

361365
# Function calling to get best-fit degree of polinomial
362366
from scipy.interpolate import interp1d
@@ -365,7 +369,7 @@ def curve_postDEM(res_name, res_directory):
365369
interpolated_elevation_values = interpolation_function(new_area_values)
366370

367371
# Function calling to get best-fit E-A-S curve values
368-
EA_curve = np.column_stack((interpolated_elevation_values, new_area_values))
372+
EA_curve = np.column_stack((interpolated_elevation_values[2:], new_area_values[2:]))
369373

370374
updated_results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]
371375
pre_area = 0
@@ -396,10 +400,10 @@ def curve_postDEM(res_name, res_directory):
396400
# Set labels and title
397401
plt.xlabel('Level (m)')
398402
plt.ylabel('Storage (mcm)')
399-
plt.title(res_name)
403+
plt.title(res_name + ' (Minimum DEM level= '+ str(round(data[0,0]))+'m)')
400404
plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
401405

402-
return min_dem
406+
return round(data[0,0])
403407

404408

405409
#=================================== Update Landsat-based E-A database with DEM-based E-A-S relationship
@@ -420,7 +424,7 @@ def one_tile(res_name, max_wl, dead_wl, res_directory):
420424
Wsa.dem_value_m[i] = closest_elev
421425
Wsa.Tot_res_volume_mcm[i] = closest_vol
422426

423-
delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] < dead_wl) | (Wsa['dem_value_m'] > max_wl+20)].index
427+
delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] > max_wl+20)].index #(Wsa['dem_value_m'] < dead_wl)
424428
Wsa = Wsa.drop(delete_id)
425429
# ========================================== EXPORT RESULTS AS A CSV FILE
426430
Wsa.to_csv('WSA_updated_' + res_name + '.csv', index=False)

‎code/InfeRes.py‎

Lines changed: 42 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -15,33 +15,38 @@
1515
from codes.CURVE import curve_preDEM
1616
from codes.CURVE import curve_postDEM
1717
from codes.MASK import mask
18-
from codes.WSA import wsa
19-
20-
18+
from codes.WSA import wsa
19+
import pandas as pd
20+
import numpy as np
21+
df = pd.read_csv('processing_res_list.csv', parse_dates=True)
2122

2223
if __name__ == "__main__":
2324

25+
2426
#====================================>> USER INPUT PARAMETERS
27+
i=0
28+
#for i in range(0, np.size(df,0)):
2529
parent_directory = "G:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
2630
os.chdir(parent_directory)
27-
res_name = "Chulabhorn"
28-
res_built_year = 1972
31+
res_name = df.Name[i]
32+
res_built_year = df.Year[i]
2933
dem_acquisition_year = 2000 #SRTM DEM (30m) acquired in Feb 2000
3034

3135
# Other Reservaoir information
3236
res_directory = parent_directory + res_name
3337
# A point within the reservoir [longitude, latitude]
34-
point = [101.636, 16.539]
38+
point = [float(value) for value in df.Point[i].split(',')]
3539
# Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
36-
boundary = [101.590, 16.617, 101.660, 16.524]
37-
max_wl = 768
40+
boundary = [float(value) for value in df.Boundary[i].split(',')]
41+
max_wl = df.Max_wl[i]
3842
os.makedirs(res_name, exist_ok=True)
3943
os.chdir(parent_directory + res_name)
4044
# Create a new folder within the working directory to download the data
4145
os.makedirs("Outputs", exist_ok=True)
4246
# Path to DEM (SouthEastAsia_DEM30m.tif), PCS: WGS1984
4347
# We have a bigger DEM file that is being used to clip the reservoir-DEM
4448
dem_file_path = "G:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif"
49+
print('Name of the reservoir: ' + res_name)
4550

4651
#------------------->> FUNCTION CALLING -1
4752
# [1]. Data pre-processing (reprojection and clipping)
@@ -61,58 +66,59 @@
6166

6267

6368
# CASE1- Reservoir built before DEM acquisition ==============================================
64-
if res_built_year < dem_acquisition_year:
69+
if res_built_year <= dem_acquisition_year:
6570
print('Name of the reservoir: ' + res_name)
6671
print('Reservoir has built before the acquisition of DEM')
6772

68-
#------------------->> FUNCTION CALLING -4
69-
# [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
70-
os.chdir(parent_directory + res_name + '/Outputs')
71-
res_minElev = curve_preDEM(res_name, point_loc, res_directory)
72-
7373
#------------------->> FUNCTION CALLING -5
74-
# [5]. Calculating the water surface area
74+
# [4.1]. Calculating the water surface area
7575
os.chdir(res_directory)
7676
wsa(res_name, res_directory)
7777

78+
#------------------->> FUNCTION CALLING -4
79+
# [5.1]. DEM-Landsat-based updated Area-Elevation-Storage curve
80+
os.chdir(parent_directory + res_name + '/Outputs')
81+
res_minElev = curve_preDEM(res_name, point_loc, res_directory)
82+
83+
7884
#------------------->> FUNCTION CALLING -6
79-
# [6]. Calculating the reservoir restorage (1 tiles)
85+
# [6.1]. Calculating the reservoir restorage (1 tiles)
8086
os.chdir(res_directory)
8187
one_tile(res_name, max_wl, res_minElev, res_directory)
82-
88+
8389
# CASE2- Reservoir built after DEM acquisition ==============================================
8490
if res_built_year > dem_acquisition_year:
8591
print('Name of the reservoir: ' + res_name)
8692
print('Reservoir has built after the acquisition of DEM')
8793

8894
#------------------->> FUNCTION CALLING -4
89-
# [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
90-
os.chdir(parent_directory + res_name + '/Outputs')
91-
res_minElev = curve_postDEM(res_name, res_directory)
92-
93-
#------------------->> FUNCTION CALLING -5
94-
# [5]. Calculating the water surface area
95+
# [4.2]. Calculating the water surface area
9596
os.chdir(res_directory)
9697
wsa(res_name, res_directory)
9798

99+
#------------------->> FUNCTION CALLING -5
100+
# [4.2]. DEM-Landsat-based updated Area-Elevation-Storage curve
101+
os.chdir(parent_directory + res_name + '/Outputs')
102+
res_minElev = curve_postDEM(res_name, max_wl, res_directory)
103+
98104
#------------------->> FUNCTION CALLING -6
99-
# [6]. Calculating the reservoir restorage (1 tiles)
105+
# [6.2]. Calculating the reservoir restorage (1 tiles)
100106
os.chdir(res_directory)
101107
one_tile(res_name, max_wl, res_minElev, res_directory)
102108

103109

104-
# ================================================ Finally moving all .png/jpg files in a seperate folder for better organisation
105-
import shutil
106-
# Create a folder to store the pictures if it doesn't exist
107-
pictures_folder = "intermediate_pictures"
108-
os.makedirs(pictures_folder, exist_ok=True)
109-
# List all files in the current directory
110-
files = os.listdir()
111-
# Move all PNG files to the 'pictures' directory
112-
for file in files:
113-
if file.lower().endswith(".png"):
114-
file_path = os.path.join(os.getcwd(), file)
115-
shutil.move(file_path, os.path.join(pictures_folder, file))
110+
# [7]. ============================ Finally moving all .png/jpg files in a seperate folder for better organisation
111+
import shutil
112+
# Create a folder to store the pictures if it doesn't exist
113+
pictures_folder = "intermediate_pictures"
114+
os.makedirs(pictures_folder, exist_ok=True)
115+
# List all files in the current directory
116+
files = os.listdir()
117+
# Move all PNG files to the 'pictures' directory
118+
for file in files:
119+
if file.lower().endswith(".png"):
120+
file_path = os.path.join(os.getcwd(), file)
121+
shutil.move(file_path, os.path.join(pictures_folder, file))
116122

117123

118124

‎code/MASK.py‎

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -101,9 +101,9 @@ def mask(res_name, max_wl, point, boundary, res_directory):
101101
print("Estimating cloud fraction...")
102102
class_count = 0
103103
cloud_threshold = 80
104-
#L8band_quality_threshold= 22280
105-
#L7band_quality_threshold= 5896
106-
#L5band_quality_threshold= 5896
104+
L8band_quality_threshold= 22280
105+
L7band_quality_threshold= 5896
106+
L5band_quality_threshold= 5896
107107
# See supplemental folder (Landsat documentation) for more information
108108
os.chdir(res_directory + "/Outputs")
109109
res_iso = gdal_array.LoadFile('res_iso.tif').astype(np.float32)
@@ -117,8 +117,8 @@ def mask(res_name, max_wl, point, boundary, res_directory):
117117
bqa = gdal_array.LoadFile(filename).astype(np.float32)
118118
water_index = "Clipped_LC08_NDWI" + filename[16:]
119119
ndwi = gdal_array.LoadFile(water_index).astype(np.float32)
120-
bqa[np.where(bqa < 22280)] = 0
121-
bqa[np.where(bqa >= 22280)] = 1
120+
bqa[np.where(bqa < L8band_quality_threshold)] = 0
121+
bqa[np.where(bqa >= L8band_quality_threshold)] = 1
122122
bqa[np.where(res_iso == 0)] = 0
123123
cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
124124
print(slno)
@@ -141,8 +141,8 @@ def mask(res_name, max_wl, point, boundary, res_directory):
141141
bqa = gdal_array.LoadFile(filename).astype(np.float32)
142142
water_index = "Clipped_LE07_NDWI"+filename[16:]
143143
ndwi = gdal_array.LoadFile(water_index).astype(np.float32)
144-
bqa[np.where(bqa < 5896)] = 0
145-
bqa[np.where(bqa >= 5896)] = 1
144+
bqa[np.where(bqa < L7band_quality_threshold)] = 0
145+
bqa[np.where(bqa >= L7band_quality_threshold)] = 1
146146
bqa[np.where(res_iso == 0)] = 0
147147
cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
148148
print(slno)
@@ -165,8 +165,8 @@ def mask(res_name, max_wl, point, boundary, res_directory):
165165
bqa = gdal_array.LoadFile(filename).astype(np.float32)
166166
water_index = "Clipped_LT05_NDWI"+filename[16:]
167167
ndwi = gdal_array.LoadFile(water_index).astype(np.float32)
168-
bqa[np.where(bqa < 5896)] = 0
169-
bqa[np.where(bqa >= 5896)] = 1
168+
bqa[np.where(bqa < L5band_quality_threshold)] = 0
169+
bqa[np.where(bqa >= L5band_quality_threshold)] = 1
170170
bqa[np.where(res_iso == 0)] = 0
171171
cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
172172
print(slno)
@@ -259,7 +259,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
259259
water_px[np.where(dem_clip <= max_wl+10)] = 1
260260
water_px[np.where(dem_clip > max_wl+10)] = 0
261261
picked_wp = pick(xp, yp, water_px)
262-
dem_mask = expand(picked_wp, 1)
262+
dem_mask = expand(picked_wp, 3)
263263
#dm_sum = np.nansum(dem_mask)
264264
output = gdal_array.SaveArray(dem_mask.astype(gdal_array.numpy.float32),
265265
"DEM_Mask.tif", format="GTiff",
@@ -356,7 +356,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
356356

357357

358358

359-
# [6] ====================== CREATE EXPANDED MASK (by 1 pixels surrounding each of water pixels)
359+
# [6] ====================== CREATE EXPANDED MASK (by 3 pixels surrounding each of water pixels)
360360
print('============ [6] CREATE EXPANDED MASK ===============')
361361
print("Creating expanded mask ...")
362362
os.chdir(res_directory + "/Outputs")
@@ -366,7 +366,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
366366
mask = sum_mask
367367
mask[np.where(sum_mask <= 1)] = 0
368368
mask[np.where(sum_mask > 1)] = 1
369-
exp_mask = expand(mask, 1)
369+
exp_mask = expand(mask, 3)
370370
output = gdal_array.SaveArray(exp_mask.astype(gdal_array.numpy.float32),
371371
"Expanded_Mask.tif",
372372
format="GTiff", prototype = res_dem_file)

0 commit comments

Comments
 (0)