diff --git a/code/CURVE.py b/code/CURVE.py index e822867..c9d1a6c 100644 --- a/code/CURVE.py +++ b/code/CURVE.py @@ -41,63 +41,43 @@ def pick(c, r, mask): # (c_number, r_number, an array of 1 amd 0) def res_iso(res_name, max_wl, point, boundary, res_directory): # ===================================================== INPUT PARAMETERS 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() - - # 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 - - bbox = [left, top, right, bottom] + 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)) - 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 - plt.figure() - plt.imshow(dem_bin, cmap='jet') - plt.colorbar() - + 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) + + #------------------ Visualization 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[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)) + #------------------ Visualization + gdal_array.SaveArray(res_iso.astype(gdal_array.numpy.float32), + "res_iso.tif", format="GTiff", + prototype = res_dem_file) #============================================================== 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 + res_dem_file = (res_name + "_DEM_UTM_CLIP.tif") + res_dem = gdal_array.LoadFile(res_dem_file).astype(np.float32) + 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), + 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() @@ -121,7 +101,7 @@ def curve(res_name, res_directory): axis=0) # saving output as a csv file - with open("Curve.csv","w", newline='') as my_csv: + with open('Curve_' + res_name + '.csv',"w", newline='') as my_csv: csvWriter = csv.writer(my_csv) csvWriter.writerows(results) diff --git a/code/CURVE_Tile.py b/code/CURVE_Tile.py index b653588..7ec9e4e 100644 --- a/code/CURVE_Tile.py +++ b/code/CURVE_Tile.py @@ -2,159 +2,239 @@ # IMPORTING LIBRARY import csv import os -import utm import numpy as np import pandas as pd from osgeo import gdal, gdal_array +import matplotlib.pyplot as plt -# HELPER FUNCTION -def pick(c, r, mask): # (c_number, r_number, an array of 1 amd 0) - filled = set() - fill = set() - fill.add((c, r)) - width = mask.shape[1]-1 - height = mask.shape[0]-1 - picked = np.zeros_like(mask, dtype=np.int8) - while fill: - x, y = fill.pop() - if y == height or x == width or x < 0 or y < 0: - continue - if mask[y][x] == 1: - picked[y][x] = 1 - filled.add((x, y)) - west = (x-1, y) - east = (x+1, y) - north = (x, y-1) - south = (x, y+1) - if west not in filled: - fill.add(west) - if east not in filled: - fill.add(east) - if north not in filled: - fill.add(north) - if south not in filled: - fill.add(south) - return picked +# # HELPER FUNCTION +# def pick(c, r, mask): # (c_number, r_number, an array of 1 amd 0) +# filled = set() +# fill = set() +# fill.add((c, r)) +# width = mask.shape[1]-1 +# height = mask.shape[0]-1 +# picked = np.zeros_like(mask, dtype=np.int8) +# while fill: +# x, y = fill.pop() +# if y == height or x == width or x < 0 or y < 0: +# continue +# if mask[y][x] == 1: +# picked[y][x] = 1 +# filled.add((x, y)) +# west = (x-1, y) +# east = (x+1, y) +# north = (x, y-1) +# south = (x, y+1) +# if west not in filled: +# fill.add(west) +# if east not in filled: +# fill.add(east) +# if north not in filled: +# fill.add(north) +# if south not in filled: +# fill.add(south) +# return picked def one_tile(res_name, max_wl, dead_wl, res_directory): #========== ESTIMATING TOTAL RESERVOIR AREA USING PREVIOUSLY GENERATED CURVE FOR THE BIGGER TILE os.chdir(res_directory + "/Outputs") - bigger = pd.read_csv("Curve.csv") - bigger_landsat_wsa = pd.read_csv("WSA.csv") + curve = pd.read_csv('Curve_' + res_name + '.csv') + landsat_wsa = pd.read_csv('WSA_' + res_name + '.csv') - Wsa= bigger_landsat_wsa + Wsa= landsat_wsa Wsa["dem_value_m"] = None Wsa["Tot_res_volume_mcm"] = None for i in range(0,len(Wsa)): - diff = np.abs(bigger.iloc[:, 1] - Wsa.Fn_area[i]) + diff = np.abs(curve.iloc[:, 1] - Wsa.Fn_area[i]) closest_index = np.argmin(diff) - closest_elev = bigger.iloc[closest_index, 0] - closest_vol = bigger.iloc[closest_index, 2] + closest_elev = curve.iloc[closest_index, 0] + closest_vol = curve.iloc[closest_index, 2] 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)].index + delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] < dead_wl) | (Wsa['dem_value_m'] > max_wl+20)].index Wsa = Wsa.drop(delete_id) # ========================================== EXPORT RESULTS AS A CSV FILE - Wsa.to_csv('WSA_Complete_' + res_name + '.csv', index=False) + Wsa.to_csv('WSA_updated_' + res_name + '.csv', index=False) print("Done") print(" ") -def two_tile(res_name, max_wl, dead_wl, point_coordinates, complete_res_boundary, dem_file_path, res_directory): - # ================================================== (Complete reservoir) - # ===================================================== INPUT PARAMETERS - bc = complete_res_boundary - pc = point_coordinates - 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) - 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 +# def two_tile(res_name, max_wl, dead_wl, point, complete_res_boundary, dem_file_path, res_directory): +# # ================================================== (Complete reservoir) +# # ===================================================== INPUT PARAMETERS +# dem = gdal.Open(dem_file_path) +# os.chdir(os.getcwd() + "/Outputs") +# res_dem_file = res_name + "DEMcomplete.tif" +# dem = gdal.Translate(res_dem_file, dem, projWin = complete_res_boundary) +# dem = None + +# #------------------ Visualization +# # image = gdal_array.LoadFile(res_dem_file) +# # plt.figure() +# # plt.imshow(image, cmap='jet') +# # plt.colorbar() + +# #============================================================ GCS to UTM +# # Input and output file paths +# ref_file = "DEM_Landsat_res_iso.tif" +# target_file = res_dem_file +# output_file = res_name + "DEMcomplete_UTM.tif" + +# # Open the reference raster (bigger dem) +# ref_ds = gdal.Open(ref_file) + +# # Get the projection and geotransform from the reference raster +# ref_proj = ref_ds.GetProjection() +# ref_transform = ref_ds.GetGeoTransform() + +# # Open the raster to be reprojected (complete dem) +# input_ds = gdal.Open(target_file) + +# # Create a new raster with the same projection as the reference raster +# output_ds = gdal.Warp(output_file, input_ds, dstSRS=ref_proj, xRes=ref_transform[1], +# yRes=abs(ref_transform[5]), +# outputBounds=[ref_transform[0], +# ref_transform[3] + ref_transform[5] * input_ds.RasterYSize, +# ref_transform[0] + ref_transform[1] * input_ds.RasterXSize, +# ref_transform[3]], +# format='GTiff') + +# # Close the input and output datasets +# input_ds = None +# output_ds = None +# ref_ds = None + +# # ===================================================== INPUT PARAMETERS +# os.chdir(res_directory + "/Outputs") +# res_dem_file = (res_name + "DEMcomplete_UTM.tif") +# dem_ds = gdal.Open(res_dem_file) +# geotransform = dem_ds.GetGeoTransform() + +# # 30m nearly equal to 0.00027777778 decimal degree +# xp = abs(round((point[0]-complete_res_boundary[0])/0.00027777778)) +# yp = abs(round((point[1]-complete_res_boundary[1])/0.00027777778)) +# dem_ds = None + +# # 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 +# #------------------ Visualization +# plt.figure() +# plt.imshow(dem_bin, cmap='jet') +# plt.colorbar() + +# 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[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)) + +#======================================================================================================= + # bc = complete_res_boundary + # pc = point_coordinates + # 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) + # 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 # CREATING E-A-S RELATIONSHOP # clipping DEM by the bounding box - print("Clipping DEM by the bounding box ...") - dem = gdal.Open(dem_file_path) + # print("Clipping DEM by the bounding box ...") + # dem = gdal.Open(dem_file_path) - # Changing path to the desired reservoir - os.chdir(os.getcwd() + "/Outputs") - res_dem_file = "Complete_" + res_name+"DEM.tif" - dem = gdal.Translate(res_dem_file, dem, projWin = bbox) - dem = None + # # Changing path to the desired reservoir + # os.chdir(os.getcwd() + "/Outputs") + # res_dem_file = "Complete_" + res_name+"DEM.tif" + # dem = gdal.Translate(res_dem_file, dem, projWin = bbox) + # dem = None - # isolating the reservoir - dem_bin = gdal_array.LoadFile(res_dem_file) - dem_bin[np.where(dem_bin > curve_ext)] = 0 - dem_bin[np.where(dem_bin > 0)] = 1 - res_iso = pick(xp, yp, dem_bin) + # # isolating the reservoir + # dem_bin = gdal_array.LoadFile(res_dem_file) + # dem_bin[np.where(dem_bin > curve_ext)] = 0 + # dem_bin[np.where(dem_bin > 0)] = 1 + # res_iso = pick(xp, yp, dem_bin) - # 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) + # # 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) - # caculating reservoir surface area and storage volume - # coresponding to each water level - 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[np.where(res_dem > i)] = 0 - water_px[np.where(water_px > 0)] = 1 - area = np.nansum(water_px)*9/10000 - storage = (area + pre_area)/2 - tot_stor += storage - pre_area = area - results = np.append(results, [[level, round(area,4), round(tot_stor,4)]], - axis=0) + # # caculating reservoir surface area and storage volume + # # coresponding to each water level + # 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[np.where(res_dem > i)] = 0 + # water_px[np.where(water_px > 0)] = 1 + # area = np.nansum(water_px)*9/10000 + # storage = (area + pre_area)/2 + # tot_stor += storage + # pre_area = area + # results = np.append(results, [[level, round(area,4), round(tot_stor,4)]], + # axis=0) - # saving output as a csv file - with open("Curve_complete_res.csv","w", newline='') as my_csv: - csvWriter = csv.writer(my_csv) - csvWriter.writerows(results) + # # saving output as a csv file + # with open("Curve_complete_res.csv","w", newline='') as my_csv: + # csvWriter = csv.writer(my_csv) + # csvWriter.writerows(results) - #========== ESTIMATING TOTAL RESERVOIR AREA USING PREVIOUSLY GENERATED CURVE FOR THE BIGGER TILE + # #========== ESTIMATING TOTAL RESERVOIR AREA USING PREVIOUSLY GENERATED CURVE FOR THE BIGGER TILE - # Import E-A curve for complete reservoir - complete = pd.read_csv("Curve_complete_res.csv") + # # Import E-A curve for complete reservoir + # complete = pd.read_csv("Curve_complete_res.csv") - # Import E-A curve for bigger portion of the reservoir - bigger = pd.read_csv("Curve.csv") + # # Import E-A curve for bigger portion of the reservoir + # bigger = pd.read_csv("Curve.csv") - # Import landsat-derived water surface area for bigger portion of the reservoir - bigger_landsat_wsa = pd.read_csv("WSA.csv") + # # Import landsat-derived water surface area for bigger portion of the reservoir + # bigger_landsat_wsa = pd.read_csv("WSA.csv") - Wsa= bigger_landsat_wsa - Wsa["dem_value_m"] = None - Wsa["Tot_res_area_skm"] = None - Wsa["Tot_res_volume_mcm"] = None - for i in range(0,len(Wsa)): - diff = np.abs(bigger.iloc[:, 1] - Wsa.Fn_area[i]) - closest_index = np.argmin(diff) - closest_elev = bigger.iloc[closest_index, 0] - Wsa.dem_value_m[i] = closest_elev + # Wsa= bigger_landsat_wsa + # Wsa["dem_value_m"] = None + # Wsa["Tot_res_area_skm"] = None + # Wsa["Tot_res_volume_mcm"] = None + # for i in range(0,len(Wsa)): + # diff = np.abs(bigger.iloc[:, 1] - Wsa.Fn_area[i]) + # closest_index = np.argmin(diff) + # closest_elev = bigger.iloc[closest_index, 0] + # Wsa.dem_value_m[i] = closest_elev - # Select the corresponding row in the complete DataFrame - selected_row = complete.iloc[np.where(complete.iloc[:, 0] == closest_elev)] - Wsa.Tot_res_area_skm[i] = selected_row.iloc[0, 1] - Wsa.Tot_res_volume_mcm[i] = selected_row.iloc[0, 2] - - delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] < dead_wl) | (Wsa['dem_value_m'] > max_wl)].index - Wsa = Wsa.drop(delete_id) - # ========================================== EXPORT RESULTS AS A CSV FILE - Wsa.to_csv('WSA_Complete_' + res_name + '.csv', index=False) - print("Done") - print(" ") + # # Select the corresponding row in the complete DataFrame + # selected_row = complete.iloc[np.where(complete.iloc[:, 0] == closest_elev)] + # Wsa.Tot_res_area_skm[i] = selected_row.iloc[0, 1] + # Wsa.Tot_res_volume_mcm[i] = selected_row.iloc[0, 2] + + # delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] < dead_wl) | (Wsa['dem_value_m'] > max_wl)].index + # Wsa = Wsa.drop(delete_id) + # # ========================================== EXPORT RESULTS AS A CSV FILE + # Wsa.to_csv('WSA_Complete_' + res_name + '.csv', index=False) + # print("Done") + # print(" ") \ No newline at end of file diff --git a/code/InfeRes.py b/code/InfeRes.py index 643855e..1ce6ba8 100644 --- a/code/InfeRes.py +++ b/code/InfeRes.py @@ -24,15 +24,13 @@ #====================================>> USER INPUT PARAMETERS parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/" os.chdir(parent_directory) - res_name = "PleiKrong" # Name of the reservoir + res_name = "Xayabouri_part2" # Name of the reservoir res_directory = parent_directory + res_name # A point within the reservoir [longitude, latitude] - point = [107.872, 14.422] + point = [101.813, 19.254] # 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 + 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 @@ -50,7 +48,7 @@ #====================================>> FUNCTION CALLING -3 # [3]. Creating mask/intermediate files - mask(res_name, yearOFcommission, max_wl, point, boundary, res_directory) + mask(res_name, max_wl, point, boundary, res_directory) #====================================>> FUNCTION CALLING -4 # [4]. DEM-Landsat-based updated Area-Elevation-Storage curve @@ -58,26 +56,14 @@ #====================================>> 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) + os.chdir(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: - 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) diff --git a/code/MASK.py b/code/MASK.py index 0d7dbbd..773bb85 100644 --- a/code/MASK.py +++ b/code/MASK.py @@ -51,27 +51,15 @@ def expand(array, n): # (an array of 1 and 0, number of additional pixels) return expand -def mask(res_name, yearOFcommission, max_wl, point, boundary, res_directory): +def mask(res_name, max_wl, point, boundary, res_directory): # [1] ======================================= INPUT PARAMETERS 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() - - # 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 - - bbox = [left, top, right, bottom] + 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)) - dem_ds = None - + yp = abs(round((point[1]-boundary[1])/0.00027777778)) # [2] =============================== DELETE >80% cloudy (over the reservoir) images @@ -231,7 +219,7 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, res_directory): print('============ [4] CREATE DEM-BASED MAX WATER EXTENT MASK ===============') print("Creating DEM-based max water extent mask ...") os.chdir(res_directory + "/Outputs") - res_dem_file = res_name + "DEM_UTM_CLIP.tif" + res_dem_file = res_name + "_DEM_UTM_CLIP.tif" dem_clip = gdal_array.LoadFile(res_dem_file).astype(np.float32) water_px = dem_clip water_px[np.where(dem_clip <= max_wl+10)] = 1 @@ -260,7 +248,7 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, res_directory): print('============ [5] CREATE LANDSAT-BASED MAX WATER EXTENT MASK ===============') print("Creating Landsat-based max water extent mask ...") os.chdir(res_directory + "/Outputs") - res_dem_file = res_name + "DEM_UTM_CLIP.tif" + res_dem_file = res_name + "_DEM_UTM_CLIP.tif" dem_clip = gdal_array.LoadFile(res_dem_file).astype(np.float32) res_iso = gdal_array.LoadFile('res_iso.tif').astype(np.float32) count = dem_clip - dem_clip @@ -325,7 +313,7 @@ def mask(res_name, yearOFcommission, max_wl, point, boundary, res_directory): 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: + with open('Landsat_Mask_' + res_name + '.csv',"w", newline='') as my_csv: csvWriter = csv.writer(my_csv) csvWriter.writerows(img_list) print("Created Landsat-based max water extent mask from "+str(img_used)+" images") diff --git a/code/Multi_Tile_Complete_WSA.py b/code/Multi_Tile_Complete_WSA.py new file mode 100644 index 0000000..bcf698e --- /dev/null +++ b/code/Multi_Tile_Complete_WSA.py @@ -0,0 +1,138 @@ +#============================================================================== +# Run this code if your reserovoir contains two or more landsat tiles. +# Run this only after succesfully processing all the individual parts (landsat tiles) of the reservoir using "InfeRes.py" +# Asign input (i.e. res_name) as "ReservoirNAME". +# Where, "ReservoirNAME" is the name of the multi-tile reservoir. +# The code will automatically search for all other parts of the reservoir and compute the combined storage +# Please save different parts of the reservoir in folder 'ReservoirNAME_tile1', 'ReservoirNAME_tile2'...and so on. +#============================================================================== + +import csv +import os +os.chdir("H:/My Drive/NUSproject/ReservoirExtraction/") +import pandas as pd +import numpy as np + +if __name__ == "__main__": + + #====================================>> USER INPUT PARAMETERS + parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/" + os.chdir(parent_directory) + multi_tile_res = 'Xayabouri' + max_wl = 285 + + print("Two tiles reservoir") + directory = os.getcwd() + filtered_files = [file for file in os.listdir(directory) if multi_tile_res in file] + data_range = pd.DataFrame() + n=0 + for filename in filtered_files: + os.chdir(parent_directory) + os.chdir(filename + '/Outputs') + curve = pd.read_csv('Curve_' + filename + '.csv') + first_column = curve.iloc[:, 0] + data_range.loc[n, 'Min'] = int(first_column.min()) + data_range.loc[n, 'Max'] = int(first_column.max()) + n += 1 + filename = None + + curve_data = [["Level (m)", "Area (skm)", "Storage (mcm)"]] + min_val= int(min(data_range.iloc[:, 0])) + max_val= int(max(data_range.iloc[:, 1])) + data = pd.DataFrame({'RefElev': range(min_val, max_val+1)}) + n=0 + for filename in filtered_files: + n += 1 + os.chdir(parent_directory) + os.chdir(filename + '/Outputs') + curve = pd.read_csv('Curve_' + filename + '.csv') + val = curve.iloc[:, 0] + data["Level"+str(n)] = 0 + data["Area"+str(n)] = 0 + data["Storage"+str(n)] = 0 + i=0 + for i in range(0,len(curve)): + #print(i) + index = np.where(data.loc[:,"RefElev"] == val.iloc[i])[0] + if (len(index)>0 and n==1): + data.iloc[index,n] = curve.iloc[i,0] + data.iloc[index,n+1] = curve.iloc[i,1] + data.iloc[index,n+2] = curve.iloc[i,2] + + if (len(index)>0 and n==2): + data.iloc[index,n+2] = curve.iloc[i,0] + data.iloc[index,n+3] = curve.iloc[i,1] + data.iloc[index,n+4] = curve.iloc[i,2] + + if (len(index)>0 and n==3): + data.iloc[index,n+4] = curve.iloc[i,0] + data.iloc[index,n+5] = curve.iloc[i,1] + data.iloc[index,n+6] = curve.iloc[i,2] + + if len(filtered_files)==2: + Lev = data['RefElev'] + Ar = data['Area1'] + data['Area2'] + Vol = data['Storage1'] + data['Storage2'] + + if len(filtered_files)==3: + Lev = data['RefElev'] + Ar = data['Area1'] + data['Area2'] + data['Area3'] + Vol = data['Storage1'] + data['Storage2'] + data['Storage3'] + + curve_data = pd.concat([Lev , Ar , Vol], axis=1, keys=["Level (m)", "Area (skm)", "Storage (mcm)"]) + curve_data['Diff'] = curve_data['Storage (mcm)'] - curve_data['Storage (mcm)'].shift(1) + pos = curve_data.index[curve_data['Diff'] <= 0].tolist() + curve_data_new = curve_data.iloc[:pos[0]] + curve_data_new = curve_data_new.drop('Diff', axis=1) + curve_data_new_list = [["Level (m)", "Area (skm)", "Storage (mcm)"]] + # Initialize curve_data_new_list with the header as the first row + curve_data_new_list = [list(curve_data_new.columns)] + # Convert DataFrame rows to a list and append to curve_data_new_list + curve_data_new_list += curve_data_new.values.tolist() + + os.chdir(parent_directory) + os.chdir(filtered_files[0] + '/Outputs') + # saving output as a csv file + with open("Curve_"+ multi_tile_res +".csv","w", newline='') as my_csv: + csvWriter = csv.writer(my_csv) + csvWriter.writerows(curve_data_new_list) + + del data, Lev, Ar, Vol, curve_data, curve_data_new + + #========== ESTIMATING TOTAL RESERVOIR AREA USING PREVIOUSLY GENERATED CURVE FOR THE BIGGER TILE + + curve = pd.read_csv('Curve_' + multi_tile_res + '.csv') + landsat_wsa = pd.read_csv('WSA_updated_' + multi_tile_res + '_tile1.csv') + dead_wl = landsat_wsa['dem_value_m'][0] + + Wsa= landsat_wsa + Wsa["Fn_area"] = None + Wsa["Tot_res_volume_mcm"] = None + for i in range(0,len(Wsa)): + index = np.where(curve.iloc[:, 0] == Wsa.dem_value_m[i]) + area = curve.iloc[index[0], 1] + volume = curve.iloc[index[0], 2] + Wsa.Fn_area[i] = area.values[0] + Wsa.Tot_res_volume_mcm[i] = volume.values[0] + + delete_id = Wsa[(Wsa['Quality'] == 0) | (Wsa['dem_value_m'] < dead_wl) | (Wsa['dem_value_m'] > max_wl+20)].index + Wsa = Wsa.drop(delete_id) + # ========================================== EXPORT RESULTS AS A CSV FILE + Wsa.to_csv('WSA_updated_' + multi_tile_res + '.csv', index=False) + print("Done") + print(" ") + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/code/PREPROCESING.py b/code/PREPROCESING.py index f694c09..2a81bfb 100644 --- a/code/PREPROCESING.py +++ b/code/PREPROCESING.py @@ -104,7 +104,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path): ndwi_subset = ndwi_ds.GetRasterBand(1).ReadAsArray(xoff, yoff, width, height) # Create the subset GeoTIFF for DEM - output_dem_path = os.path.join(res_name+"DEM_UTM_CLIP.tif") + output_dem_path = os.path.join(res_name+"_DEM_UTM_CLIP.tif") driver = gdal.GetDriverByName('GTiff') output_dem_ds = driver.Create(output_dem_path, width, height, 1, gdal.GDT_Float32) output_dem_ds.SetProjection(dem_ds.GetProjection()) @@ -113,7 +113,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path): output_dem_ds = None # Create the subset GeoTIFF for NDWI - output_ndwi_path = os.path.join(res_name+"NDWI_UTM_CLIP.tif") + output_ndwi_path = os.path.join(res_name+"_NDWI_UTM_CLIP.tif") output_ndwi_ds = driver.Create(output_ndwi_path, width, height, 1, gdal.GDT_Float32) output_ndwi_ds.SetProjection(ndwi_ds.GetProjection()) output_ndwi_ds.SetGeoTransform((xmin, ndwi_transform[1], 0, ymax, 0, ndwi_transform[5])) @@ -158,7 +158,7 @@ def preprocessing(res_name, point, boundary, parent_directory, dem_file_path): plt.savefig(res_name+'_DEM.png', dpi=600, bbox_inches='tight') #------------------ Visualization - os.remove(res_name+"NDWI_UTM_CLIP.tif") + os.remove(res_name+"_NDWI_UTM_CLIP.tif") #================================= EXTRA (UTM to GCS) =================================== # # Input and output file paths diff --git a/code/Res_bbox.xlsx b/code/Res_bbox.xlsx index f820a4e..7aa61d7 100644 Binary files a/code/Res_bbox.xlsx and b/code/Res_bbox.xlsx differ diff --git a/code/WSA.py b/code/WSA.py index 00f20d9..a9850b1 100644 --- a/code/WSA.py +++ b/code/WSA.py @@ -11,9 +11,8 @@ def wsa(res_name, res_directory): # IMPROVE NDWI-BASED LANDSAT IMAGE CLASSIFICATION results = [["Landsat", "Type", "Date", "Threshold", "R_50", "N_10", "S_zone", "Quality", "Bf_area", "Af_area", "Fn_area"]] - drtr = (res_directory + "/Outputs") + drtr = (res_directory + "/Outputs") os.chdir(res_directory + "/Clip") - #os.chdir('./Landsat_'+str(LS)) directory = os.getcwd() filtered_files = [file for file in os.listdir(directory) if "NDWI" in file] slno =0 @@ -174,7 +173,7 @@ def wsa(res_name, res_directory): # ==========================================EXPORT RESULTS AS A CSV FILE print("Exporting results as a csv file ...") os.chdir(res_directory + "/Outputs") - with open("WSA.csv","w", newline='') as my_csv: + with open('WSA_' + res_name + '.csv',"w", newline='') as my_csv: csvWriter = csv.writer(my_csv) csvWriter.writerows(results) print(" ") diff --git a/code/data_download.py b/code/data_download.py index 5f9c3f5..4bf9044 100644 --- a/code/data_download.py +++ b/code/data_download.py @@ -126,10 +126,10 @@ def clip_image_to_geometry(image): if __name__ == "__main__": # Set to the current working directory - parent_directory = "G:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/" + parent_directory = "H:/My Drive/NUSproject/ReservoirExtraction/Reservoirs/" os.chdir(parent_directory) # Name of the reservoir - res_name = "DongNaiIV" + res_name = "Xayabouri_B" 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 @@ -138,15 +138,15 @@ def clip_image_to_geometry(image): print(res_name) #====================================>> USER INPUT PARAMETERS (Landsat-8 Image Specifications) - row = 52 - path = 124 + row = 46 + path = 129 Satellite = 'L8' # A point within the reservoir [longitude, latitude] - point_coordinates = [107.7333, 11.8806] + point_coordinates = [101.8061, 19.6227] # Example coordinates [longitude, latitude] - boundary_coordinates = [107.710, 11.894, 107.875, 11.837] + boundary_coordinates = [101.754, 19.925, 102.263, 19.534] collection_id = "LANDSAT/LC08/C02/T1_TOA" - start_data = '2009-06-01' + start_data = '2018-06-01' end_date = '2022-12-31' print("-----Landsat-8 data download start-----") get_landsat_images(dataFolder, Satellite, row, path, point_coordinates, boundary_coordinates, collection_id, start_data, end_date) @@ -156,7 +156,7 @@ def clip_image_to_geometry(image): #====================================>> USER INPUT PARAMETERS (Landsat-7 Image Specifications) Satellite = 'L7' collection_id = "LANDSAT/LE07/C02/T1_TOA" - start_data = '2009-06-01' + start_data = '2018-06-01' end_date = '2022-12-31' print(res_name) print("-----Landsat-7 data download start-----") @@ -167,7 +167,7 @@ def clip_image_to_geometry(image): #====================================>> USER INPUT PARAMETERS (Landsat-5 Image Specifications) Satellite = 'L5' collection_id = "LANDSAT/LT05/C02/T1_TOA" - start_data = '2009-06-01' + start_data = '2018-06-01' end_date = '2022-12-31' print(res_name) print("-----Landsat-5 data download start-----") diff --git a/code/data_processing.py b/code/data_processing.py deleted file mode 100644 index 0488196..0000000 --- a/code/data_processing.py +++ /dev/null @@ -1,97 +0,0 @@ -############+++++++++++ 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.MASK import mask -from codes.WSA import wsa -from codes.PREPROCESING import preprocesing -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 = "Xiaowan" # Name of the reservoir - # 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] #[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.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]. Data pre-processing (reprojection and clipping) - res_minElev = preprocessing(res_name, max_wl, point, boundary, dem_file_path) - - #====================================>> FUNCTION CALLING (1) - # [1]. DEM-based Area-Elevation-Storage curve - res_minElev = curve(res_name, max_wl, point, boundary, dem_file_path) - - #====================================>> FUNCTION CALLING (2) - # [2]. Creating mask/intermediate files - res_directory = parent_directory + res_name - os.chdir(res_directory) - os.makedirs("LandsatData_Clip", exist_ok=True) - mask(res_name, yearOFcommission, max_wl, point, boundary, dem_file_path, res_directory) - - #====================================>> FUNCTION CALLING (3) - # [3]. Calculating the water surface area - res_directory = parent_directory + res_name - os.chdir(res_directory) - wsa(res_name, res_directory) - - #====================================>> FUNCTION CALLING (4) - # [4]. 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) - - - - - - - - - - - - - -