Improved footprint calculation

bin/extractSmapCyclone.py

@@ -11,11 +11,14 @@ from geoalchemy2 import Geometry
 from cyclobs_orm import SimpleTrack
 import xarray
 import pandas as pd
+import geopandas as gpd
 import numpy as np
 import math
 import glob
-
 from geo_shapely import shape360
+from scipy.spatial import ConvexHull
+from shapely.geometry import Polygon
+from shapely.wkt import dumps
 
 logging.basicConfig()
 logger = logging.getLogger(os.path.basename(__file__))
@@ -27,7 +30,7 @@ def deg180(deg360):
     convert [0,360] to [-180,180]
     """
     res = np.copy(deg360)
-    res[res > 180] = res[res > 180]-360
+    res[res > 180] = res[res > 180] - 360
     return res
 
 
@@ -98,6 +101,7 @@ def extract_start_stop_measure(dataset):
 
     return min_time, max_time
 
+
 # Set NaNs where the field "minute" and "wind" of the dataset is outside the +-time_offset
 # relative to the cyclone track point date
 def set_nan_outside_time_offset(dataset, track_data, time_offset=30):
@@ -115,6 +119,18 @@ def set_nan_outside_time_offset(dataset, track_data, time_offset=30):
 
     return dataset
 
+# Extract footprint (not taking NaN data into account)
+def compute_footprint(dataset, wind_field):
+    df = dataset.to_dataframe()
+    df = df[pd.isna(df[wind_field]) == False]
+    df = df.reset_index()
+
+    gdf = gpd.GeoDataFrame(
+        df, geometry=gpd.points_from_xy(df.lon, df.lat))
+
+    p = gdf.unary_union.convex_hull
+    return p
+
 
 def extract_write_cyclone_data(dataset, kept_track_points, filename, output_path, file_date, extract_size_km=2000):
     for sid, track_point in kept_track_points.items():
@@ -186,6 +202,9 @@ def extract_write_cyclone_data(dataset, kept_track_points, filename, output_path
         # Converting lon from 0;360 range to -180;180
         sel["lon"] = deg180(sel["lon"])
 
+        footprint_polyg = compute_footprint(sel, "wind_speed")
+        sel.attrs["footprint"] = dumps(footprint_polyg)
+
         output_filename = f"RSS_smap_wind_{track_point['date'].strftime('%Y_%m_%d_%H_%M_%S')}_{track_point['sid']}" \
                           f"_v01.0.nc"
         sel.to_netcdf(os.path.join(output_path, output_filename))
@@ -234,14 +253,6 @@ if __name__ == "__main__":
 
     engine = create_engine(args.dbd, pool_size=20, max_overflow=0)
     Session = sessionmaker(bind=engine)
-    smap_tests = [
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/205/RSS_smap_wind_daily_2020_07_23_v01.0.nc",
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/206/RSS_smap_wind_daily_2020_07_24_v01.0.nc",
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/207/RSS_smap_wind_daily_2020_07_25_v01.0.nc",
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/208/RSS_smap_wind_daily_2020_07_26_v01.0.nc",
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/209/RSS_smap_wind_daily_2020_07_27_v01.0.nc",
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/210/RSS_smap_wind_daily_2020_07_28_v01.0.nc",
-        "/home/datawork-cersat-public/provider/remss/satellite/l3/smap/smap/wind/v1.0/daily/2020/211/RSS_smap_wind_daily_2020_07_29_v01.0.nc",
-    ]
+
     for f in args.input:
         process_smap_file(Session(), f, args.output)