new commit

code/CURVE.py

@@ -38,61 +38,79 @@ 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, boundary, dem_file_path): 
+def res_iso(res_name, max_wl, point, boundary, res_directory): 
     # =====================================================  INPUT PARAMETERS
-    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 + 10 # to expand the curve  
+    os.chdir(res_directory + "/Outputs")
+    res_dem_file = (res_name + "DEM_UTM_CLIP.tif")
+    dem_ds = gdal.Open(res_dem_file)   
+    geotransform = dem_ds.GetGeoTransform()
 
-    # CREATING E-A-S RELATIONSHOP
-    # clipping DEM by the bounding box
-    print("Clipping DEM by the bounding box ...") 
-    dem = gdal.Open(dem_file_path) 
+    # 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    
 
-    # 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)   #bbox=bc
-    dem = None 
+    bbox = [left, top, right, bottom]
 
+    # 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_ds = None
+
+    # CREATING E-A-S RELATIONSHOP   
     # isolating the reservoir
     dem_bin = gdal_array.LoadFile(res_dem_file)
+    dem_bin[dem_bin == 32767] = np.nan
     #------------------ 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 > 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)
+    plt.figure()
+    plt.imshow(res_iso, cmap='jet')
+    plt.colorbar()
     output = gdal_array.SaveArray(res_iso.astype(gdal_array.numpy.float32), 
                                   "res_iso.tif", format="GTiff", 
                                   prototype = res_dem_file)
     output = None
 
     # finding the lowest DEM value in the reservoir extent
     res_dem = gdal_array.LoadFile(res_dem_file)
-    res_dem[np.where(res_iso == 0)] = 9999
-    min_dem = np.nanmin(res_dem)
+    res_dem[res_dem == 32767] = np.nan
+    res_dem[np.where(res_iso == 0)] = 9999           # 9999 is any arbitrary unrealistice value
+    min_dem = int(np.nanmin(res_dem))
 
-    # caculating reservoir surface area and storage volume 
-    # coresponding to each water level
+
+#============================================================== E-A relationship
+def curve(res_name, 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")
+    res_dem = gdal_array.LoadFile(res_dem_file)
+    res_dem[res_dem == 32767] = np.nan
+    exp_mask = gdal_array.LoadFile("Expanded_Mask.tif").astype(np.float32)
+    res_dem[np.where(exp_mask == 0)] = np.nan
+    output = gdal_array.SaveArray(res_dem.astype(gdal_array.numpy.float32), 
+                                  "DEM_Landsat_res_iso.tif", 
+                                  format="GTiff", prototype = res_dem_file)
+    output = None
+    # 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
+    res_dem_updated = ("DEM_Landsat_res_iso.tif")
+
     results = [["Level (m)", "Area (skm)", "Storage (mcm)"]]
     pre_area = 0
     tot_stor = 0 
     for i in range(min_dem, curve_ext): 
         level = i
-        water_px = gdal_array.LoadFile(res_dem_file)
-        water_px[np.where(res_iso == 0)] = 9999
+        water_px = gdal_array.LoadFile(res_dem_updated)
         water_px[np.where(res_dem > i)] = 0 
         water_px[np.where(water_px > 0)] = 1
         area = np.nansum(water_px)*9/10000
@@ -107,8 +125,7 @@ def curve(res_name, max_wl, point, boundary, dem_file_path):
         csvWriter = csv.writer(my_csv)
         csvWriter.writerows(results)
 
-# ==================== Plot the DEM-based Level-Storage curve
-
+    # ==================== Plot the DEM-based Level-Storage curve   
     data = results[1:, :]
     data = np.array(data, dtype=np.float32)
     # Extract column 2 and 3 from the array
@@ -124,6 +141,8 @@ def curve(res_name, max_wl, point, boundary, dem_file_path):
     plt.savefig(res_name+'_storageVSelevation.png', dpi=600, bbox_inches='tight')
 
     return min_dem
+
+
 
 
 

code/InfeRes.py

@@ -0,0 +1,94 @@
+############+++++++++++ 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 
+############++++++++++++++++++++++++++++++++++++++++ END ++++++++++++++++++++++++++++++++++++++++#############
+
+# IMPORTING LIBRARY
+
+import os 
+os.chdir("H:/My Drive/NUSproject/ReservoirExtraction/")
+from codes.CURVE import curve
+from codes.CURVE import res_iso
+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/"
+    os.chdir(parent_directory)
+    res_name = "PleiKrong"                        # Name of the reservoir
+    res_directory = parent_directory + res_name
+    # A point within the reservoir [longitude, latitude]
+    point = [107.872, 14.422]
+    # Upper-Left and Lower-right coordinates. Example coordinates [longitude, latitude]
+    boundary = [107.763, 14.672, 107.924, 14.392] #[107.763, 14.672, 107.924, 14.392]
+    max_wl = 575 
+    yearOFcommission = 2008
+    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)
+    # 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, yearOFcommission, 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
+    res_directory = parent_directory + res_name
+    os.chdir(res_directory)
+    wsa(res_name, res_directory)
+
+    #====================================>> FUNCTION CALLING -6
+    # [6]. Calculating the reservoir restorage (1 tiles)
+    if Number_of_tiles==1:
+        res_directory = parent_directory + res_name
+        os.chdir(res_directory)
+        print("One tile reservoir")
+        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:
+        res_directory = parent_directory + res_name
+        os.chdir(res_directory)
+        print("Two tiles reservoir")
+        # Upper-Left and Lower-right coordinates of the complete reservoir
+        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)
+
+
+
+
+
+
+
+
+
+
+
+
+
+