diff --git a/code/CURVE.py b/code/CURVE.py index cb6627a..75038ea 100644 --- a/code/CURVE.py +++ b/code/CURVE.py @@ -38,19 +38,21 @@ def pick(c, r, mask): # (c_number, r_number, an array of 1 amd 0) fill.add(south) return picked -def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_path): +def curve(res_name, max_wl, point, boundary, dem_file_path): # ===================================================== INPUT PARAMETERS - bc = boundary_coordinates - pc = point_coordinates + bc = boundary + pc = point coords = bc + pc utm_coords = np.array([utm.from_latlon(coords[i + 1], coords[i]) for i in range(0, len(coords), 2)]) # Bounding box of the reservoir [ulx, uly, lrx, lry] bbox = np.array([utm_coords[0,0], utm_coords[0,1], utm_coords[1,0], utm_coords[1,1]], dtype=np.float32) res_point = np.array([utm_coords[2,0], utm_coords[2,1]], dtype=np.float32) + # 30m equivalent distance in degree + #dist30m= 0.01745329252 xp = round(abs(res_point[0]-bbox[0])/30) yp = round(abs(res_point[1]-bbox[1])/30) # Maximum reservoir water level - curve_ext = max_wl + 20 # to expand the curve + curve_ext = max_wl + 10 # to expand the curve # CREATING E-A-S RELATIONSHOP # clipping DEM by the bounding box @@ -60,11 +62,15 @@ def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_pa # Changing path to the desired reservoir os.chdir(os.getcwd() + "/Outputs") res_dem_file = res_name+"DEM.tif" - dem = gdal.Translate(res_dem_file, dem, projWin = bbox) + dem = gdal.Translate(res_dem_file, dem, projWin = bbox) #bbox=bc dem = None # isolating the reservoir dem_bin = gdal_array.LoadFile(res_dem_file) + #------------------ Visualization + plt.figure() + plt.imshow(dem_bin, cmap='jet') + plt.colorbar() dem_bin[np.where(dem_bin > curve_ext)] = 0 dem_bin[np.where(dem_bin > 0)] = 1 res_iso = pick(xp, yp, dem_bin) @@ -114,10 +120,10 @@ def curve(res_name, max_wl, point_coordinates, boundary_coordinates, dem_file_pa # Set labels and title plt.xlabel('Level (m)') plt.ylabel('Storage (mcm)') - plt.title('Level V/s Storage curve') - - # Display the plot - plt.show() + plt.title(res_name) + plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight') + + return min_dem diff --git a/code/CURVE_Tile.py b/code/CURVE_Tile.py index 1aa6703..b653588 100644 --- a/code/CURVE_Tile.py +++ b/code/CURVE_Tile.py @@ -78,7 +78,7 @@ def two_tile(res_name, max_wl, dead_wl, point_coordinates, complete_res_boundary xp = round(abs(res_point[0]-bbox[0])/30) yp = round(abs(res_point[1]-bbox[1])/30) # Maximum reservoir water level - curve_ext = max_wl + 20 # to expand the curve + curve_ext = max_wl + 10 # to expand the curve # CREATING E-A-S RELATIONSHOP # clipping DEM by the bounding box diff --git a/code/MASK.py b/code/MASK.py index a916704..767ddd3 100644 --- a/code/MASK.py +++ b/code/MASK.py @@ -4,6 +4,7 @@ import utm import numpy as np from osgeo import gdal, gdal_array +from osgeo import osr import matplotlib.pyplot as plt #from sklearn.cluster import KMeans @@ -61,25 +62,23 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res bbox = np.array([utm_coords[0,0], utm_coords[0,1], utm_coords[1,0], utm_coords[1,1]], dtype=np.float32) res_point = np.array([utm_coords[2,0], utm_coords[2,1]], dtype=np.float32) xp = round(abs(res_point[0]-bbox[0])/30) - yp = round(abs(res_point[1]-bbox[1])/30) - - - + yp = round(abs(res_point[1]-bbox[1])/30) # [1] ============================= CLIP LANDSAT IMAGES BY THE UTM BOUNDING BOX print('============ [1] CLIP LANDSAT IMAGES BY THE UTM BOUNDING BOX ===============') print("Clipping Landsat images by the bounding box ...") clip_count = 0 - os.chdir(res_directory + "/LandsatData") + os.chdir(res_directory + "/" + res_name + "_LandsatData") directory = os.getcwd() for filename in os.listdir(directory): try: if 'BQA' in filename: year= int(filename[9:13]) if 'NDWI' in filename: year= int(filename[10:14]) - if (filename.startswith("L") and year>= yearOFcommission-2): + if (filename.startswith("L") and year>= yearOFcommission-1): + ls_img = gdal.Open(filename) print(ls_img.GetDescription()) - + output_folder = res_directory + "/LandsatData_Clip" if 'BQA' in filename: output_file = "Clipped_" + filename[:19] + '_' + res_name + filename[19:] # Output file name @@ -87,8 +86,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res output_file = "Clipped_" + filename[:20] + '_' + res_name + filename[20:] # Output file name # Create the full path for the output file - output_path = os.path.join(output_folder, output_file) - + output_path = os.path.join(output_folder, output_file) + ls_img = gdal.Translate(output_path, ls_img, projWin=bbox) ls_img = None clip_count += 1 @@ -99,7 +98,7 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res except: continue - + dem_proj=None # [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 bqa[np.where(res_iso == 0)] = 0 cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2) print(filename + " has " + str(cloud_percentage) + "% cloud coverage") - if cloud_percentage > cloud_threshold-20: + if cloud_percentage > cloud_threshold-10: print('File is removed') os.remove(filename) os.remove(water_index) @@ -265,7 +264,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res plt.figure() plt.imshow(dem_mask, cmap='jet') plt.colorbar() - plt.title('DEM_Mask.tif') + plt.title('DEM_Mask') + plt.savefig(res_name+'_DEM_Mask.png', dpi=600, bbox_inches='tight') #------------------ Visualization @@ -283,10 +283,12 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res img_list = [["Landsat", "Type", "Date"]] os.chdir(res_directory + "/LandsatData_Clip") directory = os.getcwd() - filtered_files = [file for file in os.listdir(directory) if "NDWI" in file] + filtered_filesL8 = [file for file in os.listdir(directory) if "Clipped_LC08_NDWI" in file] + filtered_filesL5 = [file for file in os.listdir(directory) if "Clipped_LT05_NDWI" in file] + filtered_files = filtered_filesL8 + filtered_filesL5 for filename in filtered_files: try: - if (filename.startswith("Clipped_LC08") or filename.startswith("Clipped_LT05")): + if (filename.startswith("Clipped_")): print(filename) ndwi = gdal_array.LoadFile(filename).astype(np.float32) ndwi[np.where(res_iso == 0)] = 0 @@ -320,7 +322,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res plt.figure() plt.imshow(count, cmap='jet') plt.colorbar() - plt.title('Count.tif') + plt.title('Count') + plt.savefig(res_name+'_Count.png', dpi=600, bbox_inches='tight') #------------------ Visualization max_we = count max_we[np.where(count < 1)] = 0 @@ -334,7 +337,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res plt.figure() plt.imshow(ls_mask, cmap='jet') plt.colorbar() - plt.title('Landsat_Mask.tif') + plt.title('Landsat_Mask') + plt.savefig(res_name+'_Landsat_Mask.png', dpi=600, bbox_inches='tight') #------------------ Visualization with open("Landsat_Mask.csv","w", newline='') as my_csv: csvWriter = csv.writer(my_csv) @@ -365,7 +369,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res plt.figure() plt.imshow(exp_mask, cmap='jet') plt.colorbar() - plt.title('Expanded_Mask.tif') + plt.title('Expanded_Mask') + plt.savefig(res_name+'_Expanded_Mask.png', dpi=600, bbox_inches='tight') #------------------ Visualization @@ -387,7 +392,8 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res plt.figure() plt.imshow(zone, cmap='jet') plt.colorbar() - plt.title('Zone_Mask.tif') + plt.title('Zone_Mask') + plt.savefig(res_name+'_Zone_Mask.png', dpi=600, bbox_inches='tight') #------------------ Visualization print("DONE!!") diff --git a/code/Res_bbox.xlsx b/code/Res_bbox.xlsx new file mode 100644 index 0000000..f820a4e Binary files /dev/null and b/code/Res_bbox.xlsx differ diff --git a/code/data_download.py b/code/data_download.py index 7cf5fce..02dca29 100644 --- a/code/data_download.py +++ b/code/data_download.py @@ -24,7 +24,7 @@ def track_task_progress(task): print("Download Completed!") # Define a function to download the satellite data from Google Earth Engine -def get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date): +def get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date): # Import the Landsat 8 TOA image collection point = ee.Geometry.Point(point_coordinates) boundary = ee.Geometry.Rectangle(boundary_coordinates) @@ -52,17 +52,22 @@ def clip_image_to_geometry(image): slno = l8images_clip.toList(l8images_clip.size()) count = 1 - for i in range(2): # range(num_images) + for i in range(num_images): # range(num_images) print(count) selected_image = ee.Image(slno.get(i)) #Change the bands as per Landsat collection [Example: B3 = Green, B5 = NIR is for Landsat8] if Satellite=='L8': ndwi = selected_image.normalizedDifference(['B3', 'B5']) - if (Satellite=='L7' or Satellite=='L5'): + #Thermal = selected_image.select('B10') + if (Satellite=='L5'): ndwi = selected_image.normalizedDifference(['B2', 'B4']) + #Thermal = selected_image.select('B6') + if (Satellite=='L7'): + ndwi = selected_image.normalizedDifference(['B2', 'B4']) + #Thermal = selected_image.select('B6_VCID_1') - BQA = selected_image.select('QA_PIXEL') + BQA = selected_image.select('QA_PIXEL') # Get the date of acquisition acquisition_date = selected_image.date() @@ -79,7 +84,7 @@ def clip_image_to_geometry(image): # Define the export parameters for a cloud-optimized GeoTIFF export_params1 = { 'image': ndwi, - 'folder': 'LandsatData', # Optional: Specify a folder in your Google Drive to save the exported image + 'folder': dataFolder, # Optional: Specify a folder in your Google Drive to save the exported image 'scale': 30, # Optional: Specify the scale/resolution of the exported image in meters 'description': collection_id[8:12] + '_NDWI_' + date_string.getInfo(), 'crs': projection['crs'], @@ -98,7 +103,7 @@ def clip_image_to_geometry(image): # Define the export parameters for a cloud-optimized GeoTIFF export_params2 = { 'image': BQA, - 'folder': 'LandsatData', # Optional: Specify a folder in your Google Drive to save the exported image + 'folder': dataFolder, # Optional: Specify a folder in your Google Drive to save the exported image 'scale': 30, # Optional: Spatial Resolution 'description': collection_id[8:12] + '_BQA_' + date_string.getInfo(), 'crs': projection['crs'], @@ -118,60 +123,69 @@ def clip_image_to_geometry(image): ################################### End of Function Definition ####################################### - if __name__ == "__main__": # Set to the current working directory parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/" os.chdir(parent_directory) # Name of the reservoir - res_name = "Juzishan" + res_name = "NamOu2" os.makedirs(res_name, exist_ok=True) os.chdir(parent_directory + res_name) # Create a new folder within the working directory to download the data - os.makedirs("LandsatData", exist_ok=True) + dataFolder = res_name + '_LandsatData' + os.makedirs(dataFolder, exist_ok=True) + print(res_name) #====================================>> USER INPUT PARAMETERS (Landsat-8 Image Specifications) - row = 43 - path = 131 + row = 46 + path = 129 Satellite = 'L8' # A point within the reservoir [longitude, latitude] - point_coordinates = [100.420, 24.66] + point_coordinates = [102.4510, 20.3926] # Example coordinates [longitude, latitude] - boundary_coordinates = [100.115, 24.79, 100.51, 24.61] + boundary_coordinates = [102.435, 20.640, 102.688, 20.388] collection_id = "LANDSAT/LC08/C02/T1_TOA" - start_data = '2001-01-01' + start_data = '2015-06-01' end_date = '2022-12-31' - get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) + print("-----Landsat-8 data download start-----") + get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) print("Congratulations...all Landsat-8 files have successfully downloaded!") + print(res_name) #====================================>> USER INPUT PARAMETERS (Landsat-7 Image Specifications) Satellite = 'L7' collection_id = "LANDSAT/LE07/C02/T1_TOA" - start_data = '2005-01-01' + start_data = '2015-06-01' end_date = '2022-12-31' - get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) + print(res_name) + print("-----Landsat-7 data download start-----") + get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) print("Congratulations...all Landsat-7 files have successfully downloaded!") + print(res_name) #====================================>> USER INPUT PARAMETERS (Landsat-5 Image Specifications) Satellite = 'L5' collection_id = "LANDSAT/LT05/C02/T1_TOA" - start_data = '2001-01-01' + start_data = '2015-06-01' end_date = '2022-12-31' - get_landsat_images(Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) + print(res_name) + print("-----Landsat-5 data download start-----") + get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) print("Congratulations...all Landsat-5 files have successfully downloaded!") + print(res_name) # ==== In case of emargency if you want to cancel the running tasks then execute the following -# # Get a list of all tasks +# Get a list of all tasks # all_tasks = ee.data.getTaskList() # # Filter out the active tasks # active_tasks = [task for task in all_tasks if task['state'] == 'RUNNING' or task['state'] == 'READY'] -# Loop through the list of active tasks and cancel each task +# #Loop through the list of active tasks and cancel each task # for task in all_tasks: # task_id = task['id'] # ee.data.cancelTask(task_id) @@ -179,3 +193,4 @@ def clip_image_to_geometry(image): + diff --git a/code/data_processing.py b/code/data_processing.py index d371bfb..14af86b 100644 --- a/code/data_processing.py +++ b/code/data_processing.py @@ -24,22 +24,28 @@ os.chdir(parent_directory) res_name = "Xiaowan" # Name of the reservoir # A point within the reservoir [longitude, latitude] - point = [100.420, 24.66] + point = [99.95, 24.745] # Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude] - boundary = [100.115, 24.79, 100.51, 24.61] - max_wl = 1020 - dead_wl = 980 - yearOFcommission = 2000 + boundary = [99.20, 25.60, 100.25, 24.65] #[107.763, 14.672, 107.924, 14.392] + max_wl = 1236 + yearOFcommission = 2010 Number_of_tiles = 1 os.makedirs(res_name, exist_ok=True) os.chdir(parent_directory + res_name) # Create a new folder within the working directory to download the data os.makedirs("Outputs", exist_ok=True) - dem_file_path = "H:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m_PCS.tif" + dem_file_path = "H:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif" + + # if boundary[2] < 102: + # dem_file_path = "G:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m_PCS_47N.tif" + # if (boundary[2] >= 102 and boundary[2] <108): + # dem_file_path = "G:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m_PCS_48N.tif" + # if (boundary[2] >= 108): + # dem_file_path = "G:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m_PCS_49N.tif" #====================================>> FUNCTION CALLING (1) # [1]. DEM-based Area-Elevation-Storage curve - curve(res_name, max_wl, point, boundary, dem_file_path) + res_minElev = curve(res_name, max_wl, point, boundary, dem_file_path) #====================================>> FUNCTION CALLING (2) # [2]. Creating mask/intermediate files @@ -60,7 +66,7 @@ res_directory = parent_directory + res_name os.chdir(res_directory) print("One tile reservoir") - one_tile(res_name, max_wl, dead_wl, res_directory) + one_tile(res_name, max_wl, res_minElev, res_directory) # Calculation of water surface area for the complete reservoir (2 tiles) and corresponding reservoir restorage if Number_of_tiles==2: @@ -68,8 +74,8 @@ os.chdir(res_directory) print("Two tiles reservoir") # Upper-Left and Lower-right coordinates of the complete reservoir - complete_res_boundary = [100.2, 23, 100.40, 22.54] - two_tile(res_name, max_wl, dead_wl, point, complete_res_boundary, dem_file_path, res_directory) + complete_res_boundary = [101.693, 20.007, 102.331, 19.240] + two_tile(res_name, max_wl, res_minElev, point, complete_res_boundary, dem_file_path, res_directory) @@ -82,7 +88,5 @@ - -