Skip to content

Commit 70f6cb5

Browse files
ssmahtossmahto
committed
commit2
1 parent 05f1661 commit 70f6cb5

File tree

6 files changed

+92
-61
lines changed

6 files changed

+92
-61
lines changed

‎code/CURVE.py‎

Lines changed: 15 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -38,19 +38,21 @@ 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_coordinates, boundary_coordinates, dem_file_path):
41+
def curve(res_name, max_wl, point, boundary, dem_file_path):
4242
# ===================================================== INPUT PARAMETERS
43-
bc = boundary_coordinates
44-
pc = point_coordinates
43+
bc = boundary
44+
pc = point
4545
coords = bc + pc
4646
utm_coords = np.array([utm.from_latlon(coords[i + 1], coords[i]) for i in range(0, len(coords), 2)])
4747
# Bounding box of the reservoir [ulx, uly, lrx, lry]
4848
bbox = np.array([utm_coords[0,0], utm_coords[0,1], utm_coords[1,0], utm_coords[1,1]], dtype=np.float32)
4949
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
5052
xp = round(abs(res_point[0]-bbox[0])/30)
5153
yp = round(abs(res_point[1]-bbox[1])/30)
5254
# Maximum reservoir water level
53-
curve_ext = max_wl + 20 # to expand the curve
55+
curve_ext = max_wl + 10 # to expand the curve
5456

5557
# CREATING E-A-S RELATIONSHOP
5658
# clipping DEM by the bounding box
@@ -60,11 +62,15 @@ def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_pa
6062
# Changing path to the desired reservoir
6163
os.chdir(os.getcwd() + "/Outputs")
6264
res_dem_file = res_name+"DEM.tif"
63-
dem = gdal.Translate(res_dem_file, dem, projWin = bbox)
65+
dem = gdal.Translate(res_dem_file, dem, projWin = bbox) #bbox=bc
6466
dem = None
6567

6668
# isolating the reservoir
6769
dem_bin = gdal_array.LoadFile(res_dem_file)
70+
#------------------ Visualization <Start>
71+
plt.figure()
72+
plt.imshow(dem_bin, cmap='jet')
73+
plt.colorbar()
6874
dem_bin[np.where(dem_bin > curve_ext)] = 0
6975
dem_bin[np.where(dem_bin > 0)] = 1
7076
res_iso = pick(xp, yp, dem_bin)
@@ -114,10 +120,10 @@ def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_pa
114120
# Set labels and title
115121
plt.xlabel('Level (m)')
116122
plt.ylabel('Storage (mcm)')
117-
plt.title('Level V/s Storage curve')
118-
119-
# Display the plot
120-
plt.show()
123+
plt.title(res_name)
124+
plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
125+
126+
return min_dem
121127

122128

123129

‎code/CURVE_Tile.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ def two_tile(res_name, max_wl, dead_wl, point_coordinates, complete_res_boundary
7878
xp = round(abs(res_point[0]-bbox[0])/30)
7979
yp = round(abs(res_point[1]-bbox[1])/30)
8080
# Maximum reservoir water level
81-
curve_ext = max_wl + 20 # to expand the curve
81+
curve_ext = max_wl + 10 # to expand the curve
8282

8383
# CREATING E-A-S RELATIONSHOP
8484
# clipping DEM by the bounding box

‎code/MASK.py‎

Lines changed: 24 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
import utm
55
import numpy as np
66
from osgeo import gdal, gdal_array
7+
from osgeo import osr
78
import matplotlib.pyplot as plt
89
#from sklearn.cluster import KMeans
910

@@ -61,34 +62,32 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
6162
bbox = np.array([utm_coords[0,0], utm_coords[0,1], utm_coords[1,0], utm_coords[1,1]], dtype=np.float32)
6263
res_point = np.array([utm_coords[2,0], utm_coords[2,1]], dtype=np.float32)
6364
xp = round(abs(res_point[0]-bbox[0])/30)
64-
yp = round(abs(res_point[1]-bbox[1])/30)
65-
66-
67-
65+
yp = round(abs(res_point[1]-bbox[1])/30)
6866

6967
# [1] ============================= CLIP LANDSAT IMAGES BY THE UTM BOUNDING BOX
7068
print('============ [1] CLIP LANDSAT IMAGES BY THE UTM BOUNDING BOX ===============')
7169
print("Clipping Landsat images by the bounding box ...")
7270
clip_count = 0
73-
os.chdir(res_directory + "/LandsatData")
71+
os.chdir(res_directory + "/" + res_name + "_LandsatData")
7472
directory = os.getcwd()
7573
for filename in os.listdir(directory):
7674
try:
7775
if 'BQA' in filename: year= int(filename[9:13])
7876
if 'NDWI' in filename: year= int(filename[10:14])
79-
if (filename.startswith("L") and year>= yearOFcommission-2):
77+
if (filename.startswith("L") and year>= yearOFcommission-1):
78+
8079
ls_img = gdal.Open(filename)
8180
print(ls_img.GetDescription())
82-
81+
8382
output_folder = res_directory + "/LandsatData_Clip"
8483
if 'BQA' in filename:
8584
output_file = "Clipped_" + filename[:19] + '_' + res_name + filename[19:] # Output file name
8685
if 'NDWI' in filename:
8786
output_file = "Clipped_" + filename[:20] + '_' + res_name + filename[20:] # Output file name
8887

8988
# Create the full path for the output file
90-
output_path = os.path.join(output_folder, output_file)
91-
89+
output_path = os.path.join(output_folder, output_file)
90+
9291
ls_img = gdal.Translate(output_path, ls_img, projWin=bbox)
9392
ls_img = None
9493
clip_count += 1
@@ -99,7 +98,7 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
9998
except:
10099
continue
101100

102-
101+
dem_proj=None
103102

104103

105104
# [2] =============================== DELETE >80% cloudy (over the reservoir) images
@@ -146,7 +145,7 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
146145
bqa[np.where(res_iso == 0)] = 0
147146
cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
148147
print(filename + " has " + str(cloud_percentage) + "% cloud coverage")
149-
if cloud_percentage > cloud_threshold-20:
148+
if cloud_percentage > cloud_threshold-10:
150149
print('File is removed')
151150
os.remove(filename)
152151
os.remove(water_index)
@@ -265,7 +264,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
265264
plt.figure()
266265
plt.imshow(dem_mask, cmap='jet')
267266
plt.colorbar()
268-
plt.title('DEM_Mask.tif')
267+
plt.title('DEM_Mask')
268+
plt.savefig(res_name+'_DEM_Mask.png', dpi=600, bbox_inches='tight')
269269
#------------------ Visualization <End>
270270

271271

@@ -283,10 +283,12 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
283283
img_list = [["Landsat", "Type", "Date"]]
284284
os.chdir(res_directory + "/LandsatData_Clip")
285285
directory = os.getcwd()
286-
filtered_files = [file for file in os.listdir(directory) if "NDWI" in file]
286+
filtered_filesL8 = [file for file in os.listdir(directory) if "Clipped_LC08_NDWI" in file]
287+
filtered_filesL5 = [file for file in os.listdir(directory) if "Clipped_LT05_NDWI" in file]
288+
filtered_files = filtered_filesL8 + filtered_filesL5
287289
for filename in filtered_files:
288290
try:
289-
if (filename.startswith("Clipped_LC08") or filename.startswith("Clipped_LT05")):
291+
if (filename.startswith("Clipped_")):
290292
print(filename)
291293
ndwi = gdal_array.LoadFile(filename).astype(np.float32)
292294
ndwi[np.where(res_iso == 0)] = 0
@@ -320,7 +322,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
320322
plt.figure()
321323
plt.imshow(count, cmap='jet')
322324
plt.colorbar()
323-
plt.title('Count.tif')
325+
plt.title('Count')
326+
plt.savefig(res_name+'_Count.png', dpi=600, bbox_inches='tight')
324327
#------------------ Visualization <End>
325328
max_we = count
326329
max_we[np.where(count < 1)] = 0
@@ -334,7 +337,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
334337
plt.figure()
335338
plt.imshow(ls_mask, cmap='jet')
336339
plt.colorbar()
337-
plt.title('Landsat_Mask.tif')
340+
plt.title('Landsat_Mask')
341+
plt.savefig(res_name+'_Landsat_Mask.png', dpi=600, bbox_inches='tight')
338342
#------------------ Visualization <End>
339343
with open("Landsat_Mask.csv","w", newline='') as my_csv:
340344
csvWriter = csv.writer(my_csv)
@@ -365,7 +369,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
365369
plt.figure()
366370
plt.imshow(exp_mask, cmap='jet')
367371
plt.colorbar()
368-
plt.title('Expanded_Mask.tif')
372+
plt.title('Expanded_Mask')
373+
plt.savefig(res_name+'_Expanded_Mask.png', dpi=600, bbox_inches='tight')
369374
#------------------ Visualization <End>
370375

371376

@@ -387,7 +392,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res
387392
plt.figure()
388393
plt.imshow(zone, cmap='jet')
389394
plt.colorbar()
390-
plt.title('Zone_Mask.tif')
395+
plt.title('Zone_Mask')
396+
plt.savefig(res_name+'_Zone_Mask.png', dpi=600, bbox_inches='tight')
391397
#------------------ Visualization <End>
392398
print("DONE!!")
393399

‎code/Res_bbox.xlsx‎

211 KB
Binary file not shown.

‎code/data_download.py‎

Lines changed: 36 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ def track_task_progress(task):
2424
print("Download Completed!")
2525

2626
# Define a function to download the satellite data from Google Earth Engine
27-
def get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date):
27+
def get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date):
2828
# Import the Landsat 8 TOA image collection
2929
point = ee.Geometry.Point(point_coordinates)
3030
boundary = ee.Geometry.Rectangle(boundary_coordinates)
@@ -52,17 +52,22 @@ def clip_image_to_geometry(image):
5252

5353
slno = l8images_clip.toList(l8images_clip.size())
5454
count = 1
55-
for i in range(2): # range(num_images)
55+
for i in range(num_images): # range(num_images)
5656
print(count)
5757
selected_image = ee.Image(slno.get(i))
5858

5959
#Change the bands as per Landsat collection [Example: B3 = Green, B5 = NIR is for Landsat8]
6060
if Satellite=='L8':
6161
ndwi = selected_image.normalizedDifference(['B3', 'B5'])
62-
if (Satellite=='L7' or Satellite=='L5'):
62+
#Thermal = selected_image.select('B10')
63+
if (Satellite=='L5'):
6364
ndwi = selected_image.normalizedDifference(['B2', 'B4'])
65+
#Thermal = selected_image.select('B6')
66+
if (Satellite=='L7'):
67+
ndwi = selected_image.normalizedDifference(['B2', 'B4'])
68+
#Thermal = selected_image.select('B6_VCID_1')
6469

65-
BQA = selected_image.select('QA_PIXEL')
70+
BQA = selected_image.select('QA_PIXEL')
6671

6772
# Get the date of acquisition
6873
acquisition_date = selected_image.date()
@@ -79,7 +84,7 @@ def clip_image_to_geometry(image):
7984
# Define the export parameters for a cloud-optimized GeoTIFF
8085
export_params1 = {
8186
'image': ndwi,
82-
'folder': 'LandsatData', # Optional: Specify a folder in your Google Drive to save the exported image
87+
'folder': dataFolder, # Optional: Specify a folder in your Google Drive to save the exported image
8388
'scale': 30, # Optional: Specify the scale/resolution of the exported image in meters
8489
'description': collection_id[8:12] + '_NDWI_' + date_string.getInfo(),
8590
'crs': projection['crs'],
@@ -98,7 +103,7 @@ def clip_image_to_geometry(image):
98103
# Define the export parameters for a cloud-optimized GeoTIFF
99104
export_params2 = {
100105
'image': BQA,
101-
'folder': 'LandsatData', # Optional: Specify a folder in your Google Drive to save the exported image
106+
'folder': dataFolder, # Optional: Specify a folder in your Google Drive to save the exported image
102107
'scale': 30, # Optional: Spatial Resolution
103108
'description': collection_id[8:12] + '_BQA_' + date_string.getInfo(),
104109
'crs': projection['crs'],
@@ -118,64 +123,74 @@ def clip_image_to_geometry(image):
118123
################################### End of Function Definition #######################################
119124

120125

121-
122126
if __name__ == "__main__":
123127

124128
# Set to the current working directory
125129
parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
126130
os.chdir(parent_directory)
127131
# Name of the reservoir
128-
res_name = "Juzishan"
132+
res_name = "NamOu2"
129133
os.makedirs(res_name, exist_ok=True)
130134
os.chdir(parent_directory + res_name)
131135
# Create a new folder within the working directory to download the data
132-
os.makedirs("LandsatData", exist_ok=True)
136+
dataFolder = res_name + '_LandsatData'
137+
os.makedirs(dataFolder, exist_ok=True)
138+
print(res_name)
133139

134140
#====================================>> USER INPUT PARAMETERS (Landsat-8 Image Specifications)
135-
row = 43
136-
path = 131
141+
row = 46
142+
path = 129
137143
Satellite = 'L8'
138144
# A point within the reservoir [longitude, latitude]
139-
point_coordinates = [100.420, 24.66]
145+
point_coordinates = [102.4510, 20.3926]
140146
# Example coordinates [longitude, latitude]
141-
boundary_coordinates = [100.115, 24.79, 100.51, 24.61]
147+
boundary_coordinates = [102.435, 20.640, 102.688, 20.388]
142148
collection_id = "LANDSAT/LC08/C02/T1_TOA"
143-
start_data = '2001-01-01'
149+
start_data = '2015-06-01'
144150
end_date = '2022-12-31'
145-
get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date)
151+
print("-----Landsat-8 data download start-----")
152+
get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date)
146153
print("Congratulations...all Landsat-8 files have successfully downloaded!")
154+
print(res_name)
147155

148156
#====================================>> USER INPUT PARAMETERS (Landsat-7 Image Specifications)
149157
Satellite = 'L7'
150158
collection_id = "LANDSAT/LE07/C02/T1_TOA"
151-
start_data = '2005-01-01'
159+
start_data = '2015-06-01'
152160
end_date = '2022-12-31'
153-
get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date)
161+
print(res_name)
162+
print("-----Landsat-7 data download start-----")
163+
get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date)
154164
print("Congratulations...all Landsat-7 files have successfully downloaded!")
165+
print(res_name)
155166

156167
#====================================>> USER INPUT PARAMETERS (Landsat-5 Image Specifications)
157168
Satellite = 'L5'
158169
collection_id = "LANDSAT/LT05/C02/T1_TOA"
159-
start_data = '2001-01-01'
170+
start_data = '2015-06-01'
160171
end_date = '2022-12-31'
161-
get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date)
172+
print(res_name)
173+
print("-----Landsat-5 data download start-----")
174+
get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date)
162175
print("Congratulations...all Landsat-5 files have successfully downloaded!")
176+
print(res_name)
163177

164178

165179

166180
# ==== In case of emargency if you want to cancel the running tasks then execute the following
167181

168-
# # Get a list of all tasks
182+
# Get a list of all tasks
169183
# all_tasks = ee.data.getTaskList()
170184

171185
# # Filter out the active tasks
172186
# active_tasks = [task for task in all_tasks if task['state'] == 'RUNNING' or task['state'] == 'READY']
173187

174-
# Loop through the list of active tasks and cancel each task
188+
# #Loop through the list of active tasks and cancel each task
175189
# for task in all_tasks:
176190
# task_id = task['id']
177191
# ee.data.cancelTask(task_id)
178192
# print(f"Cancelled task: {task_id}")
179193

180194

181195

196+

0 commit comments

Comments
 (0)