new_commit

Critical-Infrastructure-Systems-Lab · Nov 21, 2023 · 761ac08 · 761ac08
1 parent 67cdd96
commit 761ac08
Show file tree

Hide file tree

Showing 4 changed files with 164 additions and 77 deletions.
diff --git a/code/CURVE.py b/code/CURVE.py
@@ -5,7 +5,7 @@
 import os
 import utm
 import numpy as np
-from osgeo import gdal, gdal_array
+from osgeo import gdal, gdal_array, osr
 import matplotlib.pyplot as plt
 
 
@@ -38,27 +38,59 @@ def pick(c, r, mask): # (c_number, r_number, an array of 1 amd 0)
                 fill.add(south)
     return picked
 
-def res_iso(res_name, max_wl, point, boundary, res_directory): 
+
+
+def res_isolation(res_name, max_wl, point, boundary, res_directory): 
     # =====================================================  INPUT PARAMETERS
     os.chdir(res_directory + "/Outputs")
     res_dem_file = (res_name + "_DEM_UTM_CLIP.tif")
 
-    # 30m nearly equal to 0.00027777778 decimal degree
-    xp = abs(round((point[0]-boundary[0])/0.00027777778))
-    yp = abs(round((point[1]-boundary[1])/0.00027777778))     
-
-    # CREATING E-A-S RELATIONSHOP   
-    # isolating the reservoir
-    dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
-    dem_bin[dem_bin == 32767] = np.nan    
-    dem_bin[np.where(dem_bin > max_wl+10)] = 0        #to expand the reservoir extent for accounting uncertainity in max_wl
-    dem_bin[np.where(dem_bin > 0)] = 1
-    res_iso = pick(xp, yp, dem_bin)
+    try: # Converting point and boundary coordinates from CGS to UTM ===========    
+        #30m nearly equal to 0.00027777778 decimal degree
+        xp = abs(round((point[0]-boundary[0])/0.00027777778))
+        yp = abs(round((point[1]-boundary[1])/0.00027777778)) 
+
+        dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
+        dem_bin[dem_bin == 32767] = np.nan    
+        dem_bin[np.where(dem_bin > max_wl+10)] = 0        #to expand the reservoir extent for accounting uncertainity in max_wl
+        dem_bin[np.where(dem_bin > 0)] = 1
+        res_iso = pick(xp, yp, dem_bin)
+        aa=sum(sum(res_iso))
+
+        if aa == 0:    
+            dem_ds = gdal.Open(res_dem_file)   
+            geotransform = dem_ds.GetGeoTransform()   
+            # Calculate the bounding box coordinates
+            left = geotransform[0]
+            top = geotransform[3]
+            right = left + geotransform[1] * dem_ds.RasterXSize
+            bottom = top + geotransform[5] * dem_ds.RasterYSize 
+            # Bounding box of the reservoir [ulx, uly, lrx, lry]        
+            bbox = [left, top, right, bottom]
+
+            utm_coords = np.array([utm.from_latlon(point[i + 1], point[i]) for i in range(0, len(point), 2)])
+            res_point = np.array([utm_coords[0,0], utm_coords[0,1]], dtype=np.float32)
+            xp = round(abs(res_point[0]-bbox[0])/30)
+            yp = round(abs(res_point[1]-bbox[1])/30)    
+
+            dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
+            dem_bin[dem_bin == 32767] = np.nan    
+            dem_bin[np.where(dem_bin > max_wl+10)] = 0        #to expand the reservoir extent for accounting uncertainity in max_wl
+            dem_bin[np.where(dem_bin > 0)] = 1
+            res_iso = pick(xp, yp, dem_bin)
+
+    except Exception as e:
+          # Handle the exception or perform actions to handle the error gracefully
+          print(f"An error occurred: {str(e)}")
 
     #------------------ Visualization <Start>
     plt.figure()
-    plt.imshow(res_iso, cmap='jet')
-    plt.colorbar()
+    plt.imshow(res_iso, cmap='viridis')
+    plt.scatter([xp], [yp], c='r', s=10)
+    plt.imshow(dem_bin, cmap='viridis')
+    plt.scatter([xp], [yp], c='r', s=20)
+    plt.title('DEM-based reservoir isolation')
+    plt.savefig(res_name+"_DEM_res_iso.png", dpi=600, bbox_inches='tight')
     #------------------ Visualization <End>
 
     gdal_array.SaveArray(res_iso.astype(gdal_array.numpy.float32), 
@@ -81,8 +113,9 @@ def curve(res_name, res_directory):
     # plt.figure()
     # plt.imshow(res_dem, cmap='jet')
     # plt.colorbar()
+    # 
     min_dem = int(np.nanmin(res_dem))
-    curve_ext = int(np.nanmax(res_dem)) + 10              # to expand the curve
+    curve_ext = int(np.nanmax(res_dem))            
     res_dem_updated = ("DEM_Landsat_res_iso.tif")
 
     results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]

diff --git a/code/InfeRes.py b/code/InfeRes.py
@@ -1,69 +1,87 @@
 ############+++++++++++ PLEASE MAKE SURE OF THE FOLLOWING POINTS BEFORE RUNNING THE CODE +++++++++++############
 # [1]. Data are already downloaded (Satellite images and DEM)
-# [2]. DEM should be in the projected coordinate system (unit: meters)
-# [3]. Use the same coordinates that you have used in "data_download.py"
-# [4]. All the python(scripts) files are inside ".../ReservoirExtraction/codes"
-# [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 
+# [2]. Use the same coordinates that you have used in "data_download.py"
+# [3]. All the python(scripts) files are inside ".../ReservoirExtraction/codes"
 ############++++++++++++++++++++++++++++++++++++++++ END ++++++++++++++++++++++++++++++++++++++++#############
 
 # IMPORTING LIBRARY
 
 import os 
-os.chdir("H:/My Drive/NUSproject/ReservoirExtraction/")
+os.chdir("G:/My Drive/NUSproject/ReservoirExtraction/")
 from codes.CURVE import curve
-from codes.CURVE import res_iso
+from codes.CURVE import res_isolation
 from codes.MASK import mask
 from codes.WSA import wsa
 from codes.PREPROCESING import preprocessing
 from codes.CURVE_Tile import one_tile
-from codes.CURVE_Tile import two_tile
-
 
 if __name__ == "__main__":
 
     #====================================>> USER INPUT PARAMETERS 
-    parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
+    parent_directory = "G:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
     os.chdir(parent_directory)
-    res_name = "Xayabouri_part2"                        # Name of the reservoir
-    res_directory = parent_directory + res_name
-    # A point within the reservoir [longitude, latitude]
-    point = [101.813, 19.254]
-    # Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
-    boundary = [101.754, 19.534, 101.958, 19.240] #[107.763, 14.672, 107.924, 14.392]
-    max_wl = 285                            
-    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)
-    # Path to DEM (SouthEastAsia_DEM30m.tif), PCS: WGS1984
-    dem_file_path = "H:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif"
-
-    #====================================>> FUNCTION CALLING -1
-    # [1]. Data pre-processing (reprojection and clipping)
-    preprocessing(res_name, point, boundary, parent_directory, dem_file_path)
-
-    #====================================>> FUNCTION CALLING -2
-    # [2]. DEM-based reservoir isolation
-    res_iso(res_name, max_wl, point, boundary, res_directory)
-
-    #====================================>> FUNCTION CALLING -3
-    # [3]. Creating mask/intermediate files
-    mask(res_name, max_wl, point, boundary, res_directory)
-
-    #====================================>> FUNCTION CALLING -4
-    # [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
-    res_minElev = curve(res_name, res_directory)
-
-    #====================================>> FUNCTION CALLING -5
-    # [5]. Calculating the water surface area
-    os.chdir(res_directory)
-    wsa(res_name, res_directory)
-
-    #====================================>> FUNCTION CALLING -6
-    # [6]. Calculating the reservoir restorage (1 tiles)
-    os.chdir(res_directory)
-    one_tile(res_name, max_wl, res_minElev, res_directory)
-
+    res_name = "Xiaowan" 
+    res_built_year = 2010
+    dem_acquisition_year = 2000                             #SRTM DEM (30m) acquired in Feb 2000
+
+    if res_built_year > dem_acquisition_year:
+        res_directory = parent_directory + res_name
+        # A point within the reservoir [longitude, latitude]
+        point = [99.95, 24.745]
+        # Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
+        boundary = [99.20, 25.60, 100.25, 24.65] 
+        max_wl = 1236                            
+        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)
+        # Path to DEM (SouthEastAsia_DEM30m.tif), PCS: WGS1984
+        # We have a bigger DEM file that is being used to clip the reservoir-DEM
+        dem_file_path = "G:/My Drive/NUSproject/ReservoirExtraction/SEAsia_DEM/SouthEastAsia_DEM30m.tif"
+
+        #====================================>> FUNCTION CALLING -1
+        # [1]. Data pre-processing (reprojection and clipping)
+        preprocessing(res_name, point, boundary, parent_directory, dem_file_path)
+
+        #====================================>> FUNCTION CALLING -2
+        # [2]. DEM-based reservoir isolation
+        os.chdir(parent_directory + res_name + '/Outputs')
+        res_isolation(res_name, max_wl, point, boundary, res_directory)
+
+        #====================================>> FUNCTION CALLING -3
+        # [3]. Creating mask/intermediate files
+        os.chdir(parent_directory + res_name + '/Outputs')
+        mask(res_name, max_wl, point, boundary, res_directory)
+
+        #====================================>> FUNCTION CALLING -4
+        # [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
+        os.chdir(parent_directory + res_name + '/Outputs')
+        res_minElev = curve(res_name, res_directory)
+
+        #====================================>> FUNCTION CALLING -5
+        # [5]. Calculating the water surface area
+        os.chdir(res_directory)
+        wsa(res_name, res_directory)
+
+        #====================================>> FUNCTION CALLING -6
+        # [6]. Calculating the reservoir restorage (1 tiles)
+        os.chdir(res_directory)
+        one_tile(res_name, max_wl, res_minElev, res_directory)
+
+
+# Finally  moving all .png files in a seperate folder for better organisation   
+import shutil
+# Create a folder to store the pictures if it doesn't exist
+pictures_folder = "intermediate_pictures"
+os.makedirs(pictures_folder, exist_ok=True)
+# List all files in the current directory
+files = os.listdir()
+# Move all PNG files to the 'pictures' directory
+for file in files:
+    if file.lower().endswith(".png"):
+        file_path = os.path.join(os.getcwd(), file)
+        shutil.move(file_path, os.path.join(pictures_folder, file))
+
 
 
 

diff --git a/code/MASK.py b/code/MASK.py
@@ -1,6 +1,7 @@
 # IMPORTING LIBRARY
 import os
 import csv
+import utm
 import numpy as np
 from osgeo import gdal, gdal_array
 from osgeo import osr
@@ -57,10 +58,43 @@ def mask(res_name, max_wl, point, boundary, res_directory):
     os.chdir(res_directory + "/Outputs")
     res_dem_file = (res_name + "_DEM_UTM_CLIP.tif")
 
-    # 30m nearly equal to 0.00027777778 decimal degree
-    xp = abs(round((point[0]-boundary[0])/0.00027777778))
-    yp = abs(round((point[1]-boundary[1])/0.00027777778))                                          
-
+    try: # Converting point and boundary coordinates from CGS to UTM ===========    
+        #30m nearly equal to "0.00027777778" decimal degree
+        xp = abs(round((point[0]-boundary[0])/0.00027777778))
+        yp = abs(round((point[1]-boundary[1])/0.00027777778)) 
+
+        dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
+        dem_bin[dem_bin == 32767] = np.nan    
+        dem_bin[np.where(dem_bin > max_wl+10)] = 0        #to expand the reservoir extent for accounting uncertainity in max_wl
+        dem_bin[np.where(dem_bin > 0)] = 1
+        res_iso = pick(xp, yp, dem_bin)
+        aa=sum(sum(res_iso))
+
+        if aa == 0:    
+            dem_ds = gdal.Open(res_dem_file)   
+            geotransform = dem_ds.GetGeoTransform()   
+            # Calculate the bounding box coordinates
+            left = geotransform[0]
+            top = geotransform[3]
+            right = left + geotransform[1] * dem_ds.RasterXSize
+            bottom = top + geotransform[5] * dem_ds.RasterYSize 
+            # Bounding box of the reservoir [ulx, uly, lrx, lry]        
+            bbox = [left, top, right, bottom]
+
+            utm_coords = np.array([utm.from_latlon(point[i + 1], point[i]) for i in range(0, len(point), 2)])
+            res_point = np.array([utm_coords[0,0], utm_coords[0,1]], dtype=np.float32)
+            xp = round(abs(res_point[0]-bbox[0])/30)
+            yp = round(abs(res_point[1]-bbox[1])/30)    
+
+            dem_bin = gdal_array.LoadFile(res_dem_file).astype(np.float32)
+            dem_bin[dem_bin == 32767] = np.nan    
+            dem_bin[np.where(dem_bin > max_wl+10)] = 0        #to expand the reservoir extent for accounting uncertainity in max_wl
+            dem_bin[np.where(dem_bin > 0)] = 1
+            res_iso = pick(xp, yp, dem_bin)
+
+    except Exception as e:
+          # Handle the exception or perform actions to handle the error gracefully
+          print(f"An error occurred: {str(e)}")
 
     # [2] =============================== DELETE >80% cloudy (over the reservoir) images
     print('============ [2] DELETE >80% cloudy (over the reservoir) images ===============')
@@ -225,7 +259,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
     water_px[np.where(dem_clip <= max_wl+10)] = 1
     water_px[np.where(dem_clip > max_wl+10)] = 0
     picked_wp = pick(xp, yp, water_px)
-    dem_mask = expand(picked_wp, 2)
+    dem_mask = expand(picked_wp, 1)
     #dm_sum = np.nansum(dem_mask)     
     output = gdal_array.SaveArray(dem_mask.astype(gdal_array.numpy.float32), 
                                   "DEM_Mask.tif", format="GTiff", 
@@ -313,6 +347,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
     plt.title('Landsat_Mask')
     plt.savefig(res_name+'_Landsat_Mask.png', dpi=600, bbox_inches='tight')
     #------------------ Visualization <End>
+
     with open('Landsat_Mask_' + res_name + '.csv',"w", newline='') as my_csv:
         csvWriter = csv.writer(my_csv)
         csvWriter.writerows(img_list)
@@ -321,7 +356,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
 
 
 
-    # [6] ======================  CREATE EXPANDED MASK (by 2 pixels surrounding each of water pixels)
+    # [6] ======================  CREATE EXPANDED MASK (by 1 pixels surrounding each of water pixels)
     print('============ [6] CREATE EXPANDED MASK ===============')
     print("Creating expanded mask ...")
     os.chdir(res_directory +  "/Outputs")
@@ -331,7 +366,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
     mask = sum_mask
     mask[np.where(sum_mask <= 1)] = 0
     mask[np.where(sum_mask > 1)] = 1
-    exp_mask = expand(mask, 2) 
+    exp_mask = expand(mask, 1) 
     output = gdal_array.SaveArray(exp_mask.astype(gdal_array.numpy.float32), 
                                   "Expanded_Mask.tif", 
                                   format="GTiff", prototype = res_dem_file)

diff --git a/code/PREPROCESING.py b/code/PREPROCESING.py
@@ -21,7 +21,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path):
 
     # Changing path to the desired reservoir
     os.chdir(os.getcwd() + "/Outputs")
-    res_dem_file = res_name+"DEM.tif"
+    res_dem_file = res_name+"_DEM_GCS.tif"
     dem = gdal.Translate(res_dem_file, dem, projWin = boundary)
     dem = None
 
@@ -41,7 +41,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path):
     # Input and output file paths
     ref_file = (data_folder_path + '/' + first_ndwi_file)
     target_file = res_dem_file
-    output_file = res_name+"DEM_UTM.tif"     
+    output_file = res_name+"_DEM_UTM.tif"     
 
     # Open the reference raster (ndwi.tif)
     ref_ds = gdal.Open(ref_file)
@@ -151,11 +151,12 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path):
     os.chdir(parent_directory + res_name + '/Outputs')
     #------------------ Visualization <Start>
     plt.figure()
-    dem = gdal_array.LoadFile(res_name+"DEM.tif").astype(np.float32)
+    dem = gdal_array.LoadFile(res_name+"_DEM_UTM_CLIP.tif").astype(np.float32)
+    dem[dem == 32767] = np.nan
     plt.imshow(dem, cmap='jet')
     plt.colorbar()
-    plt.title('DEM')
-    plt.savefig(res_name+'_DEM.png', dpi=600, bbox_inches='tight')
+    plt.title('Clipped DEM (UTM)')
+    plt.savefig(res_name+"_DEM.png", dpi=600, bbox_inches='tight')
     #------------------ Visualization <End>
 
     os.remove(res_name+"_NDWI_UTM_CLIP.tif")