Skip to content

Commit 761ac08

Browse files
ssmahtossmahto
committed
new_commit
1 parent 67cdd96 commit 761ac08

File tree

4 files changed

+164
-77
lines changed

4 files changed

+164
-77
lines changed

‎code/CURVE.py‎

Lines changed: 49 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
import os
66
import utm
77
import numpy as np
8-
from osgeo import gdal, gdal_array
8+
from osgeo import gdal, gdal_array, osr
99
import matplotlib.pyplot as plt
1010

1111

@@ -38,27 +38,59 @@ 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 res_iso(res_name, max_wl, point, boundary, res_directory):
41+
42+
43+
def res_isolation(res_name, max_wl, point, boundary, res_directory):
4244
# ===================================================== INPUT PARAMETERS
4345
os.chdir(res_directory + "/Outputs")
4446
res_dem_file = (res_name + "_DEM_UTM_CLIP.tif")
4547

46-
# 30m nearly equal to 0.00027777778 decimal degree
47-
xp = abs(round((point[0]-boundary[0])/0.00027777778))
48-
yp = abs(round((point[1]-boundary[1])/0.00027777778))
49-
50-
# CREATING E-A-S RELATIONSHOP
51-
# isolating the reservoir
52-
dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
53-
dem_bin[dem_bin == 32767] = np.nan
54-
dem_bin[np.where(dem_bin > max_wl+10)] = 0 #to expand the reservoir extent for accounting uncertainity in max_wl
55-
dem_bin[np.where(dem_bin > 0)] = 1
56-
res_iso = pick(xp, yp, dem_bin)
48+
try: # Converting point and boundary coordinates from CGS to UTM ===========
49+
#30m nearly equal to 0.00027777778 decimal degree
50+
xp = abs(round((point[0]-boundary[0])/0.00027777778))
51+
yp = abs(round((point[1]-boundary[1])/0.00027777778))
52+
53+
dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
54+
dem_bin[dem_bin == 32767] = np.nan
55+
dem_bin[np.where(dem_bin > max_wl+10)] = 0 #to expand the reservoir extent for accounting uncertainity in max_wl
56+
dem_bin[np.where(dem_bin > 0)] = 1
57+
res_iso = pick(xp, yp, dem_bin)
58+
aa=sum(sum(res_iso))
59+
60+
if aa == 0:
61+
dem_ds = gdal.Open(res_dem_file)
62+
geotransform = dem_ds.GetGeoTransform()
63+
# Calculate the bounding box coordinates
64+
left = geotransform[0]
65+
top = geotransform[3]
66+
right = left + geotransform[1] * dem_ds.RasterXSize
67+
bottom = top + geotransform[5] * dem_ds.RasterYSize
68+
# Bounding box of the reservoir [ulx, uly, lrx, lry]
69+
bbox = [left, top, right, bottom]
70+
71+
utm_coords = np.array([utm.from_latlon(point[i + 1], point[i]) for i in range(0, len(point), 2)])
72+
res_point = np.array([utm_coords[0,0], utm_coords[0,1]], dtype=np.float32)
73+
xp = round(abs(res_point[0]-bbox[0])/30)
74+
yp = round(abs(res_point[1]-bbox[1])/30)
75+
76+
dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
77+
dem_bin[dem_bin == 32767] = np.nan
78+
dem_bin[np.where(dem_bin > max_wl+10)] = 0 #to expand the reservoir extent for accounting uncertainity in max_wl
79+
dem_bin[np.where(dem_bin > 0)] = 1
80+
res_iso = pick(xp, yp, dem_bin)
81+
82+
except Exception as e:
83+
# Handle the exception or perform actions to handle the error gracefully
84+
print(f"An error occurred: {str(e)}")
5785

5886
#------------------ Visualization <Start>
5987
plt.figure()
60-
plt.imshow(res_iso, cmap='jet')
61-
plt.colorbar()
88+
plt.imshow(res_iso, cmap='viridis')
89+
plt.scatter([xp], [yp], c='r', s=10)
90+
plt.imshow(dem_bin, cmap='viridis')
91+
plt.scatter([xp], [yp], c='r', s=20)
92+
plt.title('DEM-based reservoir isolation')
93+
plt.savefig(res_name+"_DEM_res_iso.png", dpi=600, bbox_inches='tight')
6294
#------------------ Visualization <End>
6395

6496
gdal_array.SaveArray(res_iso.astype(gdal_array.numpy.float32),
@@ -81,8 +113,9 @@ def curve(res_name, res_directory):
81113
# plt.figure()
82114
# plt.imshow(res_dem, cmap='jet')
83115
# plt.colorbar()
116+
#
84117
min_dem = int(np.nanmin(res_dem))
85-
curve_ext = int(np.nanmax(res_dem)) + 10 # to expand the curve
118+
curve_ext = int(np.nanmax(res_dem))
86119
res_dem_updated = ("DEM_Landsat_res_iso.tif")
87120

88121
results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]

‎code/InfeRes.py‎

Lines changed: 67 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -1,69 +1,87 @@
11
############+++++++++++ PLEASE MAKE SURE OF THE FOLLOWING POINTS BEFORE RUNNING THE CODE +++++++++++############
22
# [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
3+
# [2]. Use the same coordinates that you have used in "data_download.py"
4+
# [3]. All the python(scripts) files are inside ".../ReservoirExtraction/codes"
75
############++++++++++++++++++++++++++++++++++++++++ END ++++++++++++++++++++++++++++++++++++++++#############
86

97
# IMPORTING LIBRARY
108

119
import os
12-
os.chdir("H:/My Drive/NUSproject/ReservoirExtraction/")
10+
os.chdir("G:/My Drive/NUSproject/ReservoirExtraction/")
1311
from codes.CURVE import curve
14-
from codes.CURVE import res_iso
12+
from codes.CURVE import res_isolation
1513
from codes.MASK import mask
1614
from codes.WSA import wsa
1715
from codes.PREPROCESING import preprocessing
1816
from codes.CURVE_Tile import one_tile
19-
from codes.CURVE_Tile import two_tile
20-
2117

2218
if __name__ == "__main__":
2319

2420
#====================================>> USER INPUT PARAMETERS
25-
parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
21+
parent_directory = "G:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
2622
os.chdir(parent_directory)
27-
res_name = "Xayabouri_part2" # Name of the reservoir
28-
res_directory = parent_directory + res_name
29-
# A point within the reservoir [longitude, latitude]
30-
point = [101.813, 19.254]
31-
# Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
32-
boundary = [101.754, 19.534, 101.958, 19.240] #[107.763, 14.672, 107.924, 14.392]
33-
max_wl = 285
34-
os.makedirs(res_name, exist_ok=True)
35-
os.chdir(parent_directory + res_name)
36-
# Create a new folder within the working directory to download the data
37-
os.makedirs("Outputs", exist_ok=True)
38-
# Path to DEM (SouthEastAsia_DEM30m.tif), PCS: WGS1984
39-
dem_file_path = "H:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif"
40-
41-
#====================================>> FUNCTION CALLING -1
42-
# [1]. Data pre-processing (reprojection and clipping)
43-
preprocessing(res_name, point, boundary, parent_directory, dem_file_path)
44-
45-
#====================================>> FUNCTION CALLING -2
46-
# [2]. DEM-based reservoir isolation
47-
res_iso(res_name, max_wl, point, boundary, res_directory)
48-
49-
#====================================>> FUNCTION CALLING -3
50-
# [3]. Creating mask/intermediate files
51-
mask(res_name, max_wl, point, boundary, res_directory)
52-
53-
#====================================>> FUNCTION CALLING -4
54-
# [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
55-
res_minElev = curve(res_name, res_directory)
56-
57-
#====================================>> FUNCTION CALLING -5
58-
# [5]. Calculating the water surface area
59-
os.chdir(res_directory)
60-
wsa(res_name, res_directory)
61-
62-
#====================================>> FUNCTION CALLING -6
63-
# [6]. Calculating the reservoir restorage (1 tiles)
64-
os.chdir(res_directory)
65-
one_tile(res_name, max_wl, res_minElev, res_directory)
66-
23+
res_name = "Xiaowan"
24+
res_built_year = 2010
25+
dem_acquisition_year = 2000 #SRTM DEM (30m) acquired in Feb 2000
26+
27+
if res_built_year > dem_acquisition_year:
28+
res_directory = parent_directory + res_name
29+
# A point within the reservoir [longitude, latitude]
30+
point = [99.95, 24.745]
31+
# Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
32+
boundary = [99.20, 25.60, 100.25, 24.65]
33+
max_wl = 1236
34+
os.makedirs(res_name, exist_ok=True)
35+
os.chdir(parent_directory + res_name)
36+
# Create a new folder within the working directory to download the data
37+
os.makedirs("Outputs", exist_ok=True)
38+
# Path to DEM (SouthEastAsia_DEM30m.tif), PCS: WGS1984
39+
# We have a bigger DEM file that is being used to clip the reservoir-DEM
40+
dem_file_path = "G:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif"
41+
42+
#====================================>> FUNCTION CALLING -1
43+
# [1]. Data pre-processing (reprojection and clipping)
44+
preprocessing(res_name, point, boundary, parent_directory, dem_file_path)
45+
46+
#====================================>> FUNCTION CALLING -2
47+
# [2]. DEM-based reservoir isolation
48+
os.chdir(parent_directory + res_name + '/Outputs')
49+
res_isolation(res_name, max_wl, point, boundary, res_directory)
50+
51+
#====================================>> FUNCTION CALLING -3
52+
# [3]. Creating mask/intermediate files
53+
os.chdir(parent_directory + res_name + '/Outputs')
54+
mask(res_name, max_wl, point, boundary, res_directory)
55+
56+
#====================================>> FUNCTION CALLING -4
57+
# [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
58+
os.chdir(parent_directory + res_name + '/Outputs')
59+
res_minElev = curve(res_name, res_directory)
60+
61+
#====================================>> FUNCTION CALLING -5
62+
# [5]. Calculating the water surface area
63+
os.chdir(res_directory)
64+
wsa(res_name, res_directory)
65+
66+
#====================================>> FUNCTION CALLING -6
67+
# [6]. Calculating the reservoir restorage (1 tiles)
68+
os.chdir(res_directory)
69+
one_tile(res_name, max_wl, res_minElev, res_directory)
70+
71+
72+
# Finally moving all .png files in a seperate folder for better organisation
73+
import shutil
74+
# Create a folder to store the pictures if it doesn't exist
75+
pictures_folder = "intermediate_pictures"
76+
os.makedirs(pictures_folder, exist_ok=True)
77+
# List all files in the current directory
78+
files = os.listdir()
79+
# Move all PNG files to the 'pictures' directory
80+
for file in files:
81+
if file.lower().endswith(".png"):
82+
file_path = os.path.join(os.getcwd(), file)
83+
shutil.move(file_path, os.path.join(pictures_folder, file))
84+
6785

6886

6987

‎code/MASK.py‎

Lines changed: 42 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
# IMPORTING LIBRARY
22
import os
33
import csv
4+
import utm
45
import numpy as np
56
from osgeo import gdal, gdal_array
67
from osgeo import osr
@@ -57,10 +58,43 @@ def mask(res_name, max_wl, point, boundary, res_directory):
5758
os.chdir(res_directory + "/Outputs")
5859
res_dem_file = (res_name + "_DEM_UTM_CLIP.tif")
5960

60-
# 30m nearly equal to 0.00027777778 decimal degree
61-
xp = abs(round((point[0]-boundary[0])/0.00027777778))
62-
yp = abs(round((point[1]-boundary[1])/0.00027777778))
63-
61+
try: # Converting point and boundary coordinates from CGS to UTM ===========
62+
#30m nearly equal to "0.00027777778" decimal degree
63+
xp = abs(round((point[0]-boundary[0])/0.00027777778))
64+
yp = abs(round((point[1]-boundary[1])/0.00027777778))
65+
66+
dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
67+
dem_bin[dem_bin == 32767] = np.nan
68+
dem_bin[np.where(dem_bin > max_wl+10)] = 0 #to expand the reservoir extent for accounting uncertainity in max_wl
69+
dem_bin[np.where(dem_bin > 0)] = 1
70+
res_iso = pick(xp, yp, dem_bin)
71+
aa=sum(sum(res_iso))
72+
73+
if aa == 0:
74+
dem_ds = gdal.Open(res_dem_file)
75+
geotransform = dem_ds.GetGeoTransform()
76+
# Calculate the bounding box coordinates
77+
left = geotransform[0]
78+
top = geotransform[3]
79+
right = left + geotransform[1] * dem_ds.RasterXSize
80+
bottom = top + geotransform[5] * dem_ds.RasterYSize
81+
# Bounding box of the reservoir [ulx, uly, lrx, lry]
82+
bbox = [left, top, right, bottom]
83+
84+
utm_coords = np.array([utm.from_latlon(point[i + 1], point[i]) for i in range(0, len(point), 2)])
85+
res_point = np.array([utm_coords[0,0], utm_coords[0,1]], dtype=np.float32)
86+
xp = round(abs(res_point[0]-bbox[0])/30)
87+
yp = round(abs(res_point[1]-bbox[1])/30)
88+
89+
dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
90+
dem_bin[dem_bin == 32767] = np.nan
91+
dem_bin[np.where(dem_bin > max_wl+10)] = 0 #to expand the reservoir extent for accounting uncertainity in max_wl
92+
dem_bin[np.where(dem_bin > 0)] = 1
93+
res_iso = pick(xp, yp, dem_bin)
94+
95+
except Exception as e:
96+
# Handle the exception or perform actions to handle the error gracefully
97+
print(f"An error occurred: {str(e)}")
6498

6599
# [2] =============================== DELETE >80% cloudy (over the reservoir) images
66100
print('============ [2] DELETE >80% cloudy (over the reservoir) images ===============')
@@ -225,7 +259,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
225259
water_px[np.where(dem_clip <= max_wl+10)] = 1
226260
water_px[np.where(dem_clip > max_wl+10)] = 0
227261
picked_wp = pick(xp, yp, water_px)
228-
dem_mask = expand(picked_wp, 2)
262+
dem_mask = expand(picked_wp, 1)
229263
#dm_sum = np.nansum(dem_mask)
230264
output = gdal_array.SaveArray(dem_mask.astype(gdal_array.numpy.float32),
231265
"DEM_Mask.tif", format="GTiff",
@@ -313,6 +347,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
313347
plt.title('Landsat_Mask')
314348
plt.savefig(res_name+'_Landsat_Mask.png', dpi=600, bbox_inches='tight')
315349
#------------------ Visualization <End>
350+
316351
with open('Landsat_Mask_' + res_name + '.csv',"w", newline='') as my_csv:
317352
csvWriter = csv.writer(my_csv)
318353
csvWriter.writerows(img_list)
@@ -321,7 +356,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
321356

322357

323358

324-
# [6] ====================== CREATE EXPANDED MASK (by 2 pixels surrounding each of water pixels)
359+
# [6] ====================== CREATE EXPANDED MASK (by 1 pixels surrounding each of water pixels)
325360
print('============ [6] CREATE EXPANDED MASK ===============')
326361
print("Creating expanded mask ...")
327362
os.chdir(res_directory + "/Outputs")
@@ -331,7 +366,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
331366
mask = sum_mask
332367
mask[np.where(sum_mask <= 1)] = 0
333368
mask[np.where(sum_mask > 1)] = 1
334-
exp_mask = expand(mask, 2)
369+
exp_mask = expand(mask, 1)
335370
output = gdal_array.SaveArray(exp_mask.astype(gdal_array.numpy.float32),
336371
"Expanded_Mask.tif",
337372
format="GTiff", prototype = res_dem_file)

‎code/PREPROCESING.py‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path):
2121

2222
# Changing path to the desired reservoir
2323
os.chdir(os.getcwd() + "/Outputs")
24-
res_dem_file = res_name+"DEM.tif"
24+
res_dem_file = res_name+"_DEM_GCS.tif"
2525
dem = gdal.Translate(res_dem_file, dem, projWin = boundary)
2626
dem = None
2727

@@ -41,7 +41,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path):
4141
# Input and output file paths
4242
ref_file = (data_folder_path + '/' + first_ndwi_file)
4343
target_file = res_dem_file
44-
output_file = res_name+"DEM_UTM.tif"
44+
output_file = res_name+"_DEM_UTM.tif"
4545

4646
# Open the reference raster (ndwi.tif)
4747
ref_ds = gdal.Open(ref_file)
@@ -151,11 +151,12 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path):
151151
os.chdir(parent_directory + res_name + '/Outputs')
152152
#------------------ Visualization <Start>
153153
plt.figure()
154-
dem = gdal_array.LoadFile(res_name+"DEM.tif").astype(np.float32)
154+
dem = gdal_array.LoadFile(res_name+"_DEM_UTM_CLIP.tif").astype(np.float32)
155+
dem[dem == 32767] = np.nan
155156
plt.imshow(dem, cmap='jet')
156157
plt.colorbar()
157-
plt.title('DEM')
158-
plt.savefig(res_name+'_DEM.png', dpi=600, bbox_inches='tight')
158+
plt.title('Clipped DEM (UTM)')
159+
plt.savefig(res_name+"_DEM.png", dpi=600, bbox_inches='tight')
159160
#------------------ Visualization <End>
160161

161162
os.remove(res_name+"_NDWI_UTM_CLIP.tif")

0 commit comments

Comments
 (0)