commit2

Critical-Infrastructure-Systems-Lab · Nov 11, 2023 · 70f6cb5 · 70f6cb5
1 parent 05f1661
commit 70f6cb5
Show file tree

Hide file tree

Showing 6 changed files with 92 additions and 61 deletions.
diff --git 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 <Start>
+    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
@@ -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
@@ -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,34 +62,32 @@ 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
                 if 'NDWI' in filename:
                     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 <End>
 
 
@@ -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 <End>
     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 <End>
     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 <End>
 
 
@@ -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 <End>
     print("DONE!!")
 
diff --git a/code/Res_bbox.xlsx b/code/Res_bbox.xlsx
diff --git 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,64 +123,74 @@ 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)
 #     print(f"Cancelled task: {task_id}")
 
 
 
+