Skip to content

Commit 696f7f9

Browse files
ssmahtossmahto
committed
initial commit
1 parent ad06de0 commit 696f7f9

File tree

15 files changed

+463
-263
lines changed

15 files changed

+463
-263
lines changed

‎InfeRes_package/CURVE.py‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
import utm
77
import numpy as np
88
from osgeo import gdal, gdal_array
9+
import matplotlib.pyplot as plt
910

1011

1112
# HELPER FUNCTION
@@ -67,6 +68,10 @@ def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_pa
6768
dem_bin[np.where(dem_bin > curve_ext)] = 0
6869
dem_bin[np.where(dem_bin > 0)] = 1
6970
res_iso = pick(xp, yp, dem_bin)
71+
output = gdal_array.SaveArray(res_iso.astype(gdal_array.numpy.float32),
72+
"res_iso.tif", format="GTiff",
73+
prototype = res_dem_file)
74+
output = None
7075

7176
# finding the lowest DEM value in the reservoir extent
7277
res_dem = gdal_array.LoadFile(res_dem_file)
@@ -95,6 +100,24 @@ def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_pa
95100
with open("Curve.csv","w", newline='') as my_csv:
96101
csvWriter = csv.writer(my_csv)
97102
csvWriter.writerows(results)
103+
104+
# ==================== Plot the DEM-based Level-Storage curve
105+
106+
data = results[1:, :]
107+
data = np.array(data, dtype=np.float32)
108+
# Extract column 2 and 3 from the array
109+
column_1 = data[:, 0]
110+
column_3 = data[:, 2]
111+
# Create the scatter plot
112+
plt.figure()
113+
plt.scatter(column_1, column_3, s=5, c='red')
114+
# Set labels and title
115+
plt.xlabel('Level (m)')
116+
plt.ylabel('Storage (mcm)')
117+
plt.title('Level V/s Storage curve')
118+
119+
# Display the plot
120+
plt.show()
98121

99122

100123

‎InfeRes_package/MASK.py‎

Lines changed: 240 additions & 113 deletions
Large diffs are not rendered by default.

‎InfeRes_package/WSA.py‎

Lines changed: 137 additions & 134 deletions
Original file line numberDiff line numberDiff line change
@@ -19,152 +19,155 @@ def wsa(res_name, res_directory):
1919
slno =0
2020
no_zones = 50
2121
for filename in filtered_files:
22-
if filename.startswith("Clipped_"):
23-
slno = slno+1
24-
print(slno)
25-
print(filename)
26-
ndwi = gdal_array.LoadFile(filename).astype(np.float32)
27-
ndwi = np.nan_to_num(ndwi, nan = -0.5)
28-
29-
# Clip NDWI rasters by the expanded mask
30-
exp_mask = gdal_array.LoadFile(drtr+"\Expanded_Mask.tif").astype(np.float32)
31-
clip_ndwi = ndwi
32-
clip_ndwi[np.where(exp_mask == 0)] = -0.5
33-
34-
# K-means clustering clipped NDWI raters to 3 clusters
35-
# (water, wet non-water, and dry non-water) (1 more cluster for the value of -0.5)
36-
rows = len(clip_ndwi)
37-
columns = len(clip_ndwi[0])
38-
x = clip_ndwi.ravel()
39-
km = KMeans(n_clusters=4)
40-
km.fit(x.reshape(-1,1))
41-
z = km.cluster_centers_
42-
z1 = max(z)
43-
z2 = -1
44-
z3 = -1
45-
for i in range(0, 4):
46-
if z[i] < z1 and z[i] > z2:
47-
z2 = z[i]
48-
for i in range(0, 4):
49-
if z[i] < z2 and z[i] > z3:
50-
z3 = z[i]
51-
threshold = round(float((z1+z2)/2),3)
52-
print(" K-Means clustering threshold = "+str(threshold))
53-
# plt.figure(figsize=[30,15])
54-
# plt.hist(x, bins=200, range=[-0.49, 0.5], color='c')
55-
# plt.axvline(z[0], color='navy', linestyle='dashed', linewidth=2)
56-
# plt.axvline(z[1], color='navy', linestyle='dashed', linewidth=2)
57-
# plt.axvline(z[2], color='navy', linestyle='dashed', linewidth=2)
58-
# plt.axvline(z[3], color='navy', linestyle='dashed', linewidth=2)
59-
# plt.axvline((z1+z2)/2, color='red', linestyle='dashed', linewidth=2)
60-
# plt.axvline((z2+z3)/2, color='red', linestyle='dashed', linewidth=2)
61-
# plt.title(filename, fontsize=30)
62-
# plt.xlabel('NDWI', fontsize=30)
63-
# plt.xticks(fontsize=30)
64-
# plt.yticks(fontsize=30)
65-
# plt.show()
66-
labels = np.reshape(km.labels_,(-1,columns))
67-
water_label = 0
68-
for i in range(0, 4):
69-
if z[i] == max(z):
70-
water_label = i
71-
water_cluster = labels - labels
72-
water_cluster[np.where(labels == water_label)] = 1
73-
74-
# Assess image quality
75-
zone_mask = gdal_array.LoadFile(drtr+"\Zone_Mask.tif").astype(np.float32)
76-
count_zm = np.zeros(no_zones)
77-
for i in range(0, no_zones):
78-
count_zm[i] = np.count_nonzero(zone_mask == i+1)
79-
cluster_zone = zone_mask
80-
cluster_zone[np.where(water_cluster == 0)] = 0
81-
count_cl = np.zeros(no_zones)
82-
ratio = np.zeros(no_zones)
83-
N_10 = 0
84-
for i in range(0, no_zones):
85-
count_cl[i] = np.count_nonzero(cluster_zone == i+1)
86-
ratio[i] = count_cl[i]/(count_zm[i] + 1.0e-20)
87-
if ratio[i] >= 0.1:
88-
N_10 += 1
89-
print(" Ratio of zone 50 = "+str(round(ratio[no_zones-1],3)))
90-
print(" No. of zones having >=10% water pixels = "+str(int(N_10)))
91-
92-
# Improve image classification
93-
ratio_nm = ratio*100/(max(ratio) + 1.0e-20)
94-
x_axis = np.zeros(no_zones)
95-
for i in range(0, no_zones):
96-
x_axis[i] = i + 1
97-
xx = np.vstack((x_axis, ratio_nm)).T
98-
kkm = KMeans(n_clusters=2).fit(xx)
99-
llb = kkm.labels_
100-
minx0 = no_zones
101-
minx1 = no_zones
102-
for i in range(0, no_zones):
103-
if llb[i] == 0:
104-
if x_axis[i] < minx0:
105-
minx0 = x_axis[i]
106-
elif llb[i] == 1:
107-
if x_axis[i] < minx1:
108-
minx1 = x_axis[i]
109-
s_index = max(minx0, minx1)
110-
if minx0 == s_index:
111-
water_id = 0
112-
elif minx1 == s_index:
113-
water_id = 1
114-
# print(" Additional water pixels start from zone "+str(int(s_index)))
115-
# colors = ['navy' if x==water_id else 'lightblue' for x in llb]
116-
# plt.figure(figsize=[30,15])
117-
# plt.bar(x_axis, ratio, color=colors)
118-
# plt.ylim(top=1)
119-
# plt.axvline(x=s_index,color='red',linestyle='--')
120-
# plt.title(filename, fontsize=30)
121-
# plt.xlabel('Zone', fontsize=30)
122-
# plt.ylabel('Ratio', fontsize=30)
123-
# plt.xticks(fontsize=30)
124-
# plt.yticks(fontsize=30)
125-
# plt.show()
126-
127-
recall_zm = gdal_array.LoadFile(drtr+"\Zone_Mask.tif").astype(np.float32)
128-
add = recall_zm
129-
add[np.where(recall_zm < s_index)] = 0
130-
improved = added_cluster = water_cluster + add
131-
improved[np.where(added_cluster > 1)] = 1
132-
bf_area = np.count_nonzero(water_cluster == 1)*0.0009
133-
af_area = np.count_nonzero(improved == 1)*0.0009
134-
print(" Water surface area:")
135-
print(" Before improvement: "+str(round(bf_area,3))+" km2")
136-
print(" After improvement: "+str(round(af_area,3))+" km2")
137-
if bf_area == 0:
138-
fn_area = bf_area
139-
qual = 0
140-
print(" Image cannot be improved")
141-
else:
142-
if threshold < -0.5:
22+
try:
23+
if filename.startswith("Clipped_"):
24+
slno = slno+1
25+
print(slno)
26+
print(filename)
27+
ndwi = gdal_array.LoadFile(filename).astype(np.float32)
28+
ndwi = np.nan_to_num(ndwi, nan = -0.5)
29+
30+
# Clip NDWI rasters by the expanded mask
31+
exp_mask = gdal_array.LoadFile(drtr+"\Expanded_Mask.tif").astype(np.float32)
32+
clip_ndwi = ndwi
33+
clip_ndwi[np.where(exp_mask == 0)] = -0.5
34+
35+
# K-means clustering clipped NDWI raters to 3 clusters
36+
# (water, wet non-water, and dry non-water) (1 more cluster for the value of -0.5)
37+
rows = len(clip_ndwi)
38+
columns = len(clip_ndwi[0])
39+
x = clip_ndwi.ravel()
40+
km = KMeans(n_clusters=4, n_init=10)
41+
km.fit(x.reshape(-1,1))
42+
z = km.cluster_centers_
43+
z1 = max(z)
44+
z2 = -1
45+
z3 = -1
46+
for i in range(0, 4):
47+
if z[i] < z1 and z[i] > z2:
48+
z2 = z[i]
49+
for i in range(0, 4):
50+
if z[i] < z2 and z[i] > z3:
51+
z3 = z[i]
52+
threshold = round(float((z1+z2)/2),3)
53+
print(" K-Means clustering threshold = "+str(threshold))
54+
# plt.figure(figsize=[30,15])
55+
# plt.hist(x, bins=200, range=[-0.49, 0.5], color='c')
56+
# plt.axvline(z[0], color='navy', linestyle='dashed', linewidth=2)
57+
# plt.axvline(z[1], color='navy', linestyle='dashed', linewidth=2)
58+
# plt.axvline(z[2], color='navy', linestyle='dashed', linewidth=2)
59+
# plt.axvline(z[3], color='navy', linestyle='dashed', linewidth=2)
60+
# plt.axvline((z1+z2)/2, color='red', linestyle='dashed', linewidth=2)
61+
# plt.axvline((z2+z3)/2, color='red', linestyle='dashed', linewidth=2)
62+
# plt.title(filename, fontsize=30)
63+
# plt.xlabel('NDWI', fontsize=30)
64+
# plt.xticks(fontsize=30)
65+
# plt.yticks(fontsize=30)
66+
# plt.show()
67+
labels = np.reshape(km.labels_,(-1,columns))
68+
water_label = 0
69+
for i in range(0, 4):
70+
if z[i] == max(z):
71+
water_label = i
72+
water_cluster = labels - labels
73+
water_cluster[np.where(labels == water_label)] = 1
74+
75+
# Assess image quality
76+
zone_mask = gdal_array.LoadFile(drtr+"\Zone_Mask.tif").astype(np.float32)
77+
count_zm = np.zeros(no_zones)
78+
for i in range(0, no_zones):
79+
count_zm[i] = np.count_nonzero(zone_mask == i+1)
80+
cluster_zone = zone_mask
81+
cluster_zone[np.where(water_cluster == 0)] = 0
82+
count_cl = np.zeros(no_zones)
83+
ratio = np.zeros(no_zones)
84+
N_10 = 0
85+
for i in range(0, no_zones):
86+
count_cl[i] = np.count_nonzero(cluster_zone == i+1)
87+
ratio[i] = count_cl[i]/(count_zm[i] + 1.0e-20)
88+
if ratio[i] >= 0.1:
89+
N_10 += 1
90+
print(" Ratio of zone 50 = "+str(round(ratio[no_zones-1],3)))
91+
print(" No. of zones having >=10% water pixels = "+str(int(N_10)))
92+
93+
# Improve image classification
94+
ratio_nm = ratio*100/(max(ratio) + 1.0e-20)
95+
x_axis = np.zeros(no_zones)
96+
for i in range(0, no_zones):
97+
x_axis[i] = i + 1
98+
xx = np.vstack((x_axis, ratio_nm)).T
99+
kkm = KMeans(n_clusters=2).fit(xx)
100+
llb = kkm.labels_
101+
minx0 = no_zones
102+
minx1 = no_zones
103+
for i in range(0, no_zones):
104+
if llb[i] == 0:
105+
if x_axis[i] < minx0:
106+
minx0 = x_axis[i]
107+
elif llb[i] == 1:
108+
if x_axis[i] < minx1:
109+
minx1 = x_axis[i]
110+
s_index = max(minx0, minx1)
111+
if minx0 == s_index:
112+
water_id = 0
113+
elif minx1 == s_index:
114+
water_id = 1
115+
# print(" Additional water pixels start from zone "+str(int(s_index)))
116+
# colors = ['navy' if x==water_id else 'lightblue' for x in llb]
117+
# plt.figure(figsize=[30,15])
118+
# plt.bar(x_axis, ratio, color=colors)
119+
# plt.ylim(top=1)
120+
# plt.axvline(x=s_index,color='red',linestyle='--')
121+
# plt.title(filename, fontsize=30)
122+
# plt.xlabel('Zone', fontsize=30)
123+
# plt.ylabel('Ratio', fontsize=30)
124+
# plt.xticks(fontsize=30)
125+
# plt.yticks(fontsize=30)
126+
# plt.show()
127+
128+
recall_zm = gdal_array.LoadFile(drtr+"\Zone_Mask.tif").astype(np.float32)
129+
add = recall_zm
130+
add[np.where(recall_zm < s_index)] = 0
131+
improved = added_cluster = water_cluster + add
132+
improved[np.where(added_cluster > 1)] = 1
133+
bf_area = np.count_nonzero(water_cluster == 1)*0.0009
134+
af_area = np.count_nonzero(improved == 1)*0.0009
135+
print(" Water surface area:")
136+
print(" Before improvement: "+str(round(bf_area,3))+" km2")
137+
print(" After improvement: "+str(round(af_area,3))+" km2")
138+
if bf_area == 0:
143139
fn_area = bf_area
144140
qual = 0
145141
print(" Image cannot be improved")
146142
else:
147-
if ratio[49] == 0:
143+
if threshold < -0.5:
148144
fn_area = bf_area
149145
qual = 0
150146
print(" Image cannot be improved")
151147
else:
152-
if N_10 == 0:
148+
if ratio[49] == 0:
153149
fn_area = bf_area
154150
qual = 0
155151
print(" Image cannot be improved")
156152
else:
157-
fn_area = af_area
158-
qual = 1
159-
print(" Final area: "+str(round(fn_area,3))+" km2")
160-
print(" ")
161-
results = np.append(results, [[str(filename[8]), str(filename[11]),filename[18:28],
162-
round(threshold,3), round(ratio[49],3),
163-
int(N_10), int(s_index), int(qual),
164-
round(bf_area,3), round(af_area,3),
165-
round(fn_area,3)]], axis=0)
166-
continue
167-
else:
153+
if N_10 == 0:
154+
fn_area = bf_area
155+
qual = 0
156+
print(" Image cannot be improved")
157+
else:
158+
fn_area = af_area
159+
qual = 1
160+
print(" Final area: "+str(round(fn_area,3))+" km2")
161+
print(" ")
162+
results = np.append(results, [[str(filename[8]), str(filename[11]),filename[18:28],
163+
round(threshold,3), round(ratio[49],3),
164+
int(N_10), int(s_index), int(qual),
165+
round(bf_area,3), round(af_area,3),
166+
round(fn_area,3)]], axis=0)
167+
continue
168+
else:
169+
continue
170+
except:
168171
continue
169172

170173

‎InfeRes_package/data_processing.py‎

Lines changed: 13 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -8,9 +8,8 @@
88

99
# IMPORTING LIBRARY
1010

11-
import os
12-
parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/"
13-
os.chdir(parent_directory)
11+
import os
12+
os.chdir("H:/My Drive/NUSproject/ReservoirExtraction/")
1413
from codes.CURVE import curve
1514
from codes.MASK import mask
1615
from codes.WSA import wsa
@@ -20,15 +19,18 @@
2019

2120
if __name__ == "__main__":
2221

23-
#====================================>> USER INPUT PARAMETERS
22+
#====================================>> USER INPUT PARAMETERS
23+
parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
24+
os.chdir(parent_directory)
25+
res_name = "Xiaowan" # Name of the reservoir
2426
# A point within the reservoir [longitude, latitude]
25-
point = [100.373, 22.676]
27+
point = [99.95, 24.745]
2628
# Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
27-
boundary = [100.30, 23, 100.40, 22.54]
28-
res_name = "Nuozhadu" # Name of the reservoir
29-
max_wl = 812
30-
dead_wl = 750
31-
Number_of_tiles = 2
29+
boundary = [99.20, 25.60, 100.25, 24.65]
30+
max_wl = 1236
31+
dead_wl = 1162
32+
yearOFcommission = 2010
33+
Number_of_tiles = 1
3234
os.makedirs(res_name, exist_ok=True)
3335
os.chdir(parent_directory + res_name)
3436
# Create a new folder within the working directory to download the data
@@ -44,7 +46,7 @@
4446
res_directory = parent_directory + res_name
4547
os.chdir(res_directory)
4648
os.makedirs("LandsatData_Clip", exist_ok=True)
47-
mask(res_name, max_wl, point, boundary, dem_file_path, res_directory)
49+
mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res_directory)
4850

4951
#====================================>> FUNCTION CALLING (3)
5052
# [3]. Calculating the water surface area

‎InfeRes_package/images/Count.png‎

69.2 KB
Loading

‎InfeRes_package/images/DEMbased_Level_Storage_Curve.png‎

11.8 KB
Loading

‎InfeRes_package/images/DemMASK.png‎

36.5 KB
Loading

‎InfeRes_package/images/ExpandedMASK.png‎

36.5 KB
Loading

‎InfeRes_package/images/LandsatMASK.png‎

36.8 KB
Loading

‎InfeRes_package/images/XiaowanDEM.png‎

214 KB
Loading

0 commit comments

Comments
 (0)