new_commit

code/CURVE.py

@@ -59,7 +59,7 @@ def best_fit_degree (xdata, ydata):
     from sklearn.metrics import r2_score
 
     # Maximum degree for polynomial regression
-    max_degree = 5
+    max_degree = 4
 
     # Lists to store R-squared values and RMSEs for each degree
     r_squared_values = []
@@ -154,8 +154,8 @@ def curve_extrapolate(xdata, ydata, bf_degree):
     coefficients = np.polyfit(xdata, ydata, degree)
     p = np.poly1d(coefficients)
 
-    # Generate area values approaching zero for extrapolation
-    extrapolated_xdata_values = np.linspace(0, max(xdata), 1000)  
+    # Generate area values approaching 30% of the largest area for extrapolation
+    extrapolated_xdata_values = np.linspace(0.2*min(xdata), max(xdata), 1000)  
     #extrapolated_area_values = np.linspace(max(area), 1.5*max(area), 1000)
     # Use the polynomial function to extrapolate elevation values for small area values
     extrapolated_ydata_values = p(extrapolated_xdata_values)
@@ -267,8 +267,8 @@ def curve_preDEM(res_name, point_loc, res_directory):
     # Extract column 2 and 3 from the array    
     data = results[1:, :]
     data = np.array(data, dtype=np.float32)
-    area = data[:, 1]    # area 
-    elevation = data[:, 0]    # elevation (ydata)
+    area = data[2:, 1]    # area 
+    elevation = data[2:, 0]    # elevation (ydata)
 
     # Function calling to get best-fit degree of polinomial
     bf_degree = best_fit_degree(area, elevation)
@@ -308,15 +308,15 @@ def curve_preDEM(res_name, point_loc, res_directory):
     # Set labels and title
     plt.xlabel('Level (m)')
     plt.ylabel('Storage (mcm)')
-    plt.title(res_name)
+    plt.title(res_name + ' (Minimum DEM level= '+ str(round(data[0,0]))+'m)')
     plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
 
     return round(data[0,0])  
 
 
 
 #============================================================== E-A relationship
-def curve_postDEM(res_name, res_directory): 
+def curve_postDEM(res_name, max_wl, res_directory): 
     # caculating reservoir surface area and storage volume coresponding to each water level
     os.chdir(res_directory + "/Outputs")                    
     res_dem_file = (res_name + "_DEM_UTM_CLIP.tif")
@@ -334,7 +334,7 @@ def curve_postDEM(res_name, res_directory):
     # plt.colorbar()
     # 
     min_dem = int(np.nanmin(res_dem))
-    curve_ext = int(np.nanmax(res_dem))            
+    curve_ext = max_wl+20            
     res_dem_updated = ("DEM_Landsat_res_iso.tif")
 
     results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]
@@ -355,8 +355,12 @@ def curve_postDEM(res_name, res_directory):
     # Extract column 2 and 3 from the array    
     data = results[1:, :]
     data = np.array(data, dtype=np.float32)
-    area = data[:, 1]    # area 
-    elevation = data[:, 0]    # elevation (ydata)
+    idx = np.where(data[:, 1]==min(data[:, 1]))
+    data = data[np.size(idx,1)-1:]
+    idx = np.where((data[:, 0] > 0) & (data[:, 1] > 0) & (data[:, 2] > 0))[0]
+    data = data[idx[0]:]
+    area = data[2:, 1]    # area 
+    elevation = data[2:, 0]    # elevation (ydata)
 
     # Function calling to get best-fit degree of polinomial
     from scipy.interpolate import interp1d
@@ -365,7 +369,7 @@ def curve_postDEM(res_name, res_directory):
     interpolated_elevation_values = interpolation_function(new_area_values)
 
     # Function calling to get best-fit E-A-S curve values
-    EA_curve = np.column_stack((interpolated_elevation_values, new_area_values))  
+    EA_curve = np.column_stack((interpolated_elevation_values[2:], new_area_values[2:]))  
 
     updated_results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]
     pre_area = 0
@@ -396,10 +400,10 @@ def curve_postDEM(res_name, res_directory):
     # Set labels and title
     plt.xlabel('Level (m)')
     plt.ylabel('Storage (mcm)')
-    plt.title(res_name)
+    plt.title(res_name + ' (Minimum DEM level= '+ str(round(data[0,0]))+'m)')
     plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
 
-    return min_dem
+    return round(data[0,0])
 
 
 #=================================== Update Landsat-based E-A database with DEM-based E-A-S relationship
@@ -420,7 +424,7 @@ def one_tile(res_name, max_wl, dead_wl, res_directory):
         Wsa.dem_value_m[i] = closest_elev
         Wsa.Tot_res_volume_mcm[i] = closest_vol
 
-    delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] < dead_wl) | (Wsa['dem_value_m'] > max_wl+20)].index
+    delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] > max_wl+20)].index  #(Wsa['dem_value_m'] < dead_wl)
     Wsa = Wsa.drop(delete_id)
     # ========================================== EXPORT RESULTS AS A CSV FILE
     Wsa.to_csv('WSA_updated_' + res_name + '.csv', index=False)

code/InfeRes.py

@@ -15,33 +15,38 @@
 from codes.CURVE import curve_preDEM
 from codes.CURVE import curve_postDEM
 from codes.MASK import mask
-from codes.WSA import wsa
-
-
+from codes.WSA import wsa 
+import pandas as pd
+import numpy as np
+df = pd.read_csv('processing_res_list.csv', parse_dates=True)
 
 if __name__ == "__main__":
 
+
     #====================================>> USER INPUT PARAMETERS 
+    i=0
+    #for i in range(0, np.size(df,0)):
     parent_directory = "G:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/"
     os.chdir(parent_directory)
-    res_name = "Chulabhorn" 
-    res_built_year = 1972
+    res_name = df.Name[i] 
+    res_built_year = df.Year[i]
     dem_acquisition_year = 2000                             #SRTM DEM (30m) acquired in Feb 2000
 
     # Other Reservaoir information
     res_directory = parent_directory + res_name
     # A point within the reservoir [longitude, latitude]
-    point = [101.636, 16.539]
+    point = [float(value) for value in df.Point[i].split(',')]
     # Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
-    boundary = [101.590, 16.617, 101.660, 16.524] 
-    max_wl = 768                            
+    boundary = [float(value) for value in df.Boundary[i].split(',')] 
+    max_wl = df.Max_wl[i]                            
     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"
+    print('Name of the reservoir: ' + res_name)
 
     #------------------->> FUNCTION CALLING -1
     # [1]. Data pre-processing (reprojection and clipping)
@@ -61,58 +66,59 @@
 
 
     # CASE1- Reservoir built before DEM acquisition ==============================================
-    if res_built_year < dem_acquisition_year:      
+    if res_built_year <= dem_acquisition_year:      
         print('Name of the reservoir: ' + res_name)
         print('Reservoir has built before the acquisition of DEM')
 
-        #------------------->> FUNCTION CALLING -4
-        # [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
-        os.chdir(parent_directory + res_name + '/Outputs')
-        res_minElev = curve_preDEM(res_name, point_loc, res_directory)
-
         #------------------->> FUNCTION CALLING -5
-        # [5]. Calculating the water surface area
+        # [4.1]. Calculating the water surface area
         os.chdir(res_directory)
         wsa(res_name, res_directory)
 
+        #------------------->> FUNCTION CALLING -4
+        # [5.1]. DEM-Landsat-based updated Area-Elevation-Storage curve
+        os.chdir(parent_directory + res_name + '/Outputs')
+        res_minElev = curve_preDEM(res_name, point_loc, res_directory)
+
+
         #------------------->> FUNCTION CALLING -6
-        # [6]. Calculating the reservoir restorage (1 tiles)
+        # [6.1]. Calculating the reservoir restorage (1 tiles)
         os.chdir(res_directory)
         one_tile(res_name, max_wl, res_minElev, res_directory)
-
+    
     # CASE2- Reservoir built after DEM acquisition ==============================================
     if res_built_year > dem_acquisition_year:        
         print('Name of the reservoir: ' + res_name)
         print('Reservoir has built after the acquisition of DEM')
 
         #------------------->> FUNCTION CALLING -4
-        # [4]. DEM-Landsat-based updated Area-Elevation-Storage curve
-        os.chdir(parent_directory + res_name + '/Outputs')
-        res_minElev = curve_postDEM(res_name, res_directory)
-
-        #------------------->> FUNCTION CALLING -5
-        # [5]. Calculating the water surface area
+        # [4.2]. Calculating the water surface area
         os.chdir(res_directory)
         wsa(res_name, res_directory)
 
+        #------------------->> FUNCTION CALLING -5
+        # [4.2]. DEM-Landsat-based updated Area-Elevation-Storage curve
+        os.chdir(parent_directory + res_name + '/Outputs')
+        res_minElev = curve_postDEM(res_name, max_wl, res_directory)
+
         #------------------->> FUNCTION CALLING -6
-        # [6]. Calculating the reservoir restorage (1 tiles)
+        # [6.2]. Calculating the reservoir restorage (1 tiles)
         os.chdir(res_directory)
         one_tile(res_name, max_wl, res_minElev, res_directory)
 
 
-# ================================================ Finally  moving all .png/jpg 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))
+    # [7]. ============================ Finally  moving all .png/jpg 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))
 
 
 

code/MASK.py

@@ -101,9 +101,9 @@ def mask(res_name, max_wl, point, boundary, res_directory):
     print("Estimating cloud fraction...")
     class_count = 0 
     cloud_threshold = 80
-    #L8band_quality_threshold= 22280     
-    #L7band_quality_threshold= 5896
-    #L5band_quality_threshold= 5896
+    L8band_quality_threshold= 22280     
+    L7band_quality_threshold= 5896
+    L5band_quality_threshold= 5896
     # See supplemental folder (Landsat documentation) for more information
     os.chdir(res_directory + "/Outputs") 
     res_iso = gdal_array.LoadFile('res_iso.tif').astype(np.float32)
@@ -117,8 +117,8 @@ def mask(res_name, max_wl, point, boundary, res_directory):
                 bqa = gdal_array.LoadFile(filename).astype(np.float32)
                 water_index = "Clipped_LC08_NDWI" + filename[16:]
                 ndwi = gdal_array.LoadFile(water_index).astype(np.float32)
-                bqa[np.where(bqa < 22280)] = 0
-                bqa[np.where(bqa >= 22280)] = 1
+                bqa[np.where(bqa < L8band_quality_threshold)] = 0
+                bqa[np.where(bqa >= L8band_quality_threshold)] = 1
                 bqa[np.where(res_iso == 0)] = 0
                 cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
                 print(slno)
@@ -141,8 +141,8 @@ def mask(res_name, max_wl, point, boundary, res_directory):
                 bqa = gdal_array.LoadFile(filename).astype(np.float32)
                 water_index = "Clipped_LE07_NDWI"+filename[16:]
                 ndwi = gdal_array.LoadFile(water_index).astype(np.float32)
-                bqa[np.where(bqa < 5896)] = 0
-                bqa[np.where(bqa >= 5896)] = 1
+                bqa[np.where(bqa < L7band_quality_threshold)] = 0
+                bqa[np.where(bqa >= L7band_quality_threshold)] = 1
                 bqa[np.where(res_iso == 0)] = 0
                 cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
                 print(slno)
@@ -165,8 +165,8 @@ def mask(res_name, max_wl, point, boundary, res_directory):
                 bqa = gdal_array.LoadFile(filename).astype(np.float32)
                 water_index = "Clipped_LT05_NDWI"+filename[16:]
                 ndwi = gdal_array.LoadFile(water_index).astype(np.float32)
-                bqa[np.where(bqa < 5896)] = 0
-                bqa[np.where(bqa >= 5896)] = 1
+                bqa[np.where(bqa < L5band_quality_threshold)] = 0
+                bqa[np.where(bqa >= L5band_quality_threshold)] = 1
                 bqa[np.where(res_iso == 0)] = 0
                 cloud_percentage = round(np.sum(bqa)/np.sum(res_iso)*100,2)
                 print(slno)
@@ -259,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, 1)
+    dem_mask = expand(picked_wp, 3)
     #dm_sum = np.nansum(dem_mask)     
     output = gdal_array.SaveArray(dem_mask.astype(gdal_array.numpy.float32), 
                                   "DEM_Mask.tif", format="GTiff", 
@@ -356,7 +356,7 @@ def mask(res_name, max_wl, point, boundary, res_directory):
 
 
 
-    # [6] ======================  CREATE EXPANDED MASK (by 1 pixels surrounding each of water pixels)
+    # [6] ======================  CREATE EXPANDED MASK (by 3 pixels surrounding each of water pixels)
     print('============ [6] CREATE EXPANDED MASK ===============')
     print("Creating expanded mask ...")
     os.chdir(res_directory +  "/Outputs")
@@ -366,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, 1) 
+    exp_mask = expand(mask, 3) 
     output = gdal_array.SaveArray(exp_mask.astype(gdal_array.numpy.float32), 
                                   "Expanded_Mask.tif", 
                                   format="GTiff", prototype = res_dem_file)