write to disk

GloFAS/GloFAS_prep/station_to_glofas.py

@@ -0,0 +1,34 @@
+from comparison.pointMatching import find_closestcorresponding_point
+import pandas as pd 
+import xarray as xr 
+import GloFAS.GloFAS_prep.configuration as cfg
+
+# Load data
+stations_ups_csv = f'{cfg.DataDir}/Stations_ups_area_DNH.csv'
+upstream_area_nc = f"{cfg.DataDir}/uparea_glofas_v4_0.nc"
+
+stations_upsarea_df = pd.read_csv(stations_ups_csv)
+
+glofas_ups_area_ds = xr.open_dataset(upstream_area_nc, engine="h5netcdf")
+glofas_ups_area_da = glofas_ups_area_ds['uparea']  # Access the upstream area variable
+
+# A wrapper function for apply
+def find_closest_wrapper(row):
+    return find_closestcorresponding_point(
+        point_x=row['Lon'],
+        point_y=row['Lat'],
+        ups_area_point_m=row['Catchment area (km2)'] * 1e6,  # Convert from km² to m²
+        glofas_ups_area_xr=glofas_ups_area_da,
+        radius_m=6000
+    )
+
+# Apply the function row-wise
+results = stations_upsarea_df.apply(find_closest_wrapper, axis=1)
+
+# Extract results into separate columns
+stations_upsarea_df['Glofas_Point_X'], stations_upsarea_df['Glofas_Point_Y'], stations_upsarea_df['Area_Diff'] = zip(*results)
+
+
+# Print or save the updated DataFrame
+print(stations_upsarea_df.head())
+stations_upsarea_df.to_csv(f'{cfg.DataDir}/stations_Glofas_coordinates.csv')

comparison/pointMatching.py

@@ -34,19 +34,6 @@ def findclosestpoint(point_x, point_y, target_gdf):
     dist, idx = tree.query([(point_x, point_y)], k=1)  # Find closest point
     return target_gdf.iloc[idx[0]]  # Return row of the closest point
 
-def find_closestcorresponding_point (point_x, point_y, ups_area_point_m,  glofas_ups_area_xr, radius_m = 5000): 
-    '''within the radius_m radius, from the point_x and point_y, find the cell in glofas_ups_area_xr for which the value corresponds most to ups_area_point
-    point_x: coordinate of station x 
-    point_y: coordinate of station y 
-    ups_area_point: upstream area (also called catchment area) for that specific station in m^3
-    model_ups_area_xr: xr.DataArray entailing the upstream areas in the model m^3
-    radius_m = radius around the point for which to look for the most corresponding upstream area value
-    '''
-    import numpy as np
-import xarray as xr
-from geopy.distance import geodesic
-
-
 
 def find_closestcorresponding_point(point_x, point_y, ups_area_point_m, glofas_ups_area_xr, radius_m=5000):
     """
@@ -67,8 +54,8 @@ def find_closestcorresponding_point(point_x, point_y, ups_area_point_m, glofas_u
     """
 
     # Extract coordinates and values from the DataArray
-    lats = glofas_ups_area_xr["lat"].values
-    lons = glofas_ups_area_xr["lon"].values
+    lats = glofas_ups_area_xr["latitude"].values
+    lons = glofas_ups_area_xr["longitude"].values
     areas = glofas_ups_area_xr.values
 
     # Create a meshgrid of the coordinates