Merge pull request #46 from gregory-halverson/main

v2.11.0 robust input dataset

Author: gregory-halverson

FLiESANN/ECOv002-cal-val-FLiESANN-inputs.csv

@@ -6,10 +6,11 @@
 from .prepare_FLiESANN_inputs import prepare_FLiESANN_inputs
 from .run_FLiESANN_inference import run_FLiESANN_inference
 from .process_FLiESANN import FLiESANN
-from .generate_FLiESANN_inputs_table import generate_FLiES_inputs_table
+from .generate_FLiESANN_inputs_table_deprecated import generate_FLiES_inputs_table
 from .process_FLiESANN_table import process_FLiESANN_table
 from .ECOv002_static_tower_FLiESANN_inputs import load_ECOv002_static_tower_FLiESANN_inputs
 from .ECOv002_calval_FLiESANN_inputs import load_ECOv002_calval_FLiESANN_inputs
 from .ECOv002_calval_FLiESANN_outputs import load_ECOv002_calval_FLiESANN_outputs
 from .verify import verify
 from .retrieve_FLiESANN_GEOS5FP_inputs import retrieve_FLiESANN_GEOS5FP_inputs
+from .generate_FLiESANN_inputs_table import generate_FLiESANN_inputs_table

FLiESANN/ECOv002-cal-val-FLiESANN-outputs.csv

@@ -0,0 +1,65 @@
+import numpy as np
+import pandas as pd
+import shapely
+import rasters as rt
+
+
+def filter_dataframe_to_location_time_pairs(df, geometry, time_UTC):
+    """Filter DataFrame or 2D array returned from GEOS5FP to match original location-time pairs"""
+    
+    # Handle DataFrame
+    if isinstance(df, pd.DataFrame):
+        if len(df) == len(geometry.geoms):
+            return df
+            
+        # Extract first column (data values) from DataFrame
+        data_array = df.iloc[:, 0].values.astype(np.float32)
+    # Handle 2D numpy array
+    elif isinstance(df, np.ndarray) and len(df.shape) == 2:
+        if df.shape[0] == (len(geometry.geoms) if hasattr(geometry, 'geoms') else 0):
+            return df
+        
+        # Extract first column if multiple columns
+        if df.shape[1] > 1:
+            data_array = df[:, 0].astype(np.float32)
+        else:
+            data_array = df.flatten().astype(np.float32)
+    # Handle 1D array or scalar
+    else:
+        return df
+    
+    # Get coordinates of input geometry points
+    if isinstance(geometry, (shapely.geometry.MultiPoint, rt.MultiPoint)):
+        input_coords = [(geom.x, geom.y) for geom in geometry.geoms]
+    else:
+        return df
+    
+    # Convert time_UTC to array if it's a single value
+    if not hasattr(time_UTC, '__len__'):
+        time_UTC_array = [time_UTC] * len(input_coords)
+    else:
+        time_UTC_array = time_UTC
+    
+    # For each input location-time pair, find the matching DataFrame row
+    # GEOS5FP processes unique times and returns data for all locations at each time
+    # We need to select only the rows that match our specific location-time pairs
+    
+    # Create a mapping of (lat, lon, time_index) -> row index
+    # The DataFrame contains all locations for each unique time
+    unique_times = sorted(set(pd.to_datetime(time_UTC_array)))
+    time_to_index = {t: i for i, t in enumerate(unique_times)}
+    
+    selected_rows = []
+    for i, (coord, time_val) in enumerate(zip(input_coords, time_UTC_array)):
+        time_val = pd.to_datetime(time_val)
+        time_idx = time_to_index[time_val]
+        # Row index = time_idx * num_locations + location_idx
+        row_idx = time_idx * len(input_coords) + i
+        if row_idx < len(data_array):
+            selected_rows.append(row_idx)
+    
+    if len(selected_rows) == len(input_coords):
+        return data_array[selected_rows].astype(np.float32)
+    else:
+        # Fallback: return all rows if filtering fails
+        return data_array.astype(np.float32)

FLiESANN/FLiESANN.py

@@ -1,162 +1,139 @@
-from dateutil import parser
+import logging
 
 import numpy as np
 import pandas as pd
-
+import rasters as rt
+from dateutil import parser
+from pandas import DataFrame
+from rasters import MultiPoint, WGS84
+from shapely.geometry import Point
 from GEOS5FP import GEOS5FP
-from sentinel_tiles import sentinel_tiles
-from sun_angles import calculate_SZA_from_datetime
-from koppengeiger import load_koppen_geiger
-
-import logging
+from NASADEM import NASADEMConnection
+from .retrieve_FLiESANN_inputs import retrieve_FLiESANN_inputs
 
 logger = logging.getLogger(__name__)
 
-def generate_FLiES_inputs_table(
-        FLiES_inputs_from_calval_df: pd.DataFrame,
-        GEOS5FP_connection: GEOS5FP = None) -> pd.DataFrame:
-    """
-    FLiES_inputs_from_claval_df:
-        Pandas DataFrame containing the columns: tower, lat, lon, time_UTC, albedo, elevation_km
-    return:
-        Pandas DataFrame containing the columns: tower, lat, lon, time_UTC, doy, albedo, elevation_km, AOT, COT, vapor_gccm, ozone_cm, SZA, KG
+def generate_FLiESANN_inputs_table(
+        input_df: DataFrame,
+        GEOS5FP_connection: GEOS5FP = None,
+        NASADEM_connection: NASADEMConnection = None) -> DataFrame:
     """
-    if GEOS5FP_connection is None:
-        GEOS5FP_connection = GEOS5FP()
-
-    # output_rows = []
-    FLiES_inputs_df = FLiES_inputs_from_calval_df.copy()
-
-    doy = []
-    AOT = []
-    COT = []
-    vapor_gccm = []
-    ozone_cm = []
-    SWin = []
-    Tmin_K = []
-    SZA = []
-    KG = []
-
-    for i, input_row in FLiES_inputs_from_calval_df.iterrows():
-        tower = input_row.tower
-        lat = input_row.lat
-        lon = input_row.lon
-        time_UTC = input_row.time_UTC
-        albedo = input_row.albedo
-        elevation_km = input_row.elevation_km
-        logger.info(f"collecting FLiES inputs for tower {tower} lat {lat} lon {lon} time {time_UTC} UTC")
-        time_UTC = parser.parse(str(time_UTC))
-        doy.append(time_UTC.timetuple().tm_yday)
-        date_UTC = time_UTC.date()
-        tile = sentinel_tiles.toMGRS(lat, lon)[:5]
-
-        try:
-            tile_grid = sentinel_tiles.grid(tile=tile, cell_size=70)
-        except Exception as e:
-            logger.error(e)
-            logger.warning(f"unable to process tile {tile}")
-            AOT.append(np.nan)
-            COT.append(np.nan)
-            vapor_gccm.append(np.nan)
-            ozone_cm.append(np.nan)
-            SZA.append(np.nan)
-            KG.append(np.nan)
-            continue
-
-        rows, cols = tile_grid.shape
-        row, col = tile_grid.index_point(rt.Point(lon, lat))
-        geometry = tile_grid[max(0, row - 1):min(row + 2, rows - 1),
-                             max(0, col - 1):min(col + 2, cols - 1)]
-
-        if not "AOT" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("retrieving GEOS-5 FP aerosol optical thickness raster")
-                AOT.append(np.nanmedian(GEOS5FP_connection.AOT(time_UTC=time_UTC, geometry=geometry)))
-            except Exception as e:
-                AOT.append(np.nan)
-                logger.exception(e)
-
-        if not "COT" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("generating GEOS-5 FP cloud optical thickness raster")
-                COT.append(np.nanmedian(GEOS5FP_connection.COT(time_UTC=time_UTC, geometry=geometry)))
-            except Exception as e:
-                COT.append(np.nan)
-                logger.exception(e)
-        
-        if not "vapor_gccm" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("generating GEOS5-FP water vapor raster in grams per square centimeter")
-                vapor_gccm.append(np.nanmedian(GEOS5FP_connection.vapor_gccm(time_UTC=time_UTC, geometry=geometry)))
-            except Exception as e:
-                vapor_gccm.append(np.nan)
-                logger.exception(e)
-
-        if not "ozone_cm" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("generating GEOS5-FP ozone raster in grams per square centimeter")
-                ozone_cm.append(np.nanmedian(GEOS5FP_connection.ozone_cm(time_UTC=time_UTC, geometry=geometry)))
-            except Exception as e:
-                ozone_cm.append(np.nan)
-                logger.exception(e)
-        
-        if not "SWin" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("generating GEOS5-FP incoming solar radiation raster in watts per square meter")
-                SWin.append(np.nanmedian(GEOS5FP_connection.SWin(time_UTC=time_UTC, geometry=geometry)))
-            except Exception as e:
-                SWin.append(np.nan)
-                logger.exception(e)
-
-        if not "Tmin_K" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("generating GEOS5-FP minimum temperature in Kelvin")
-                Tmin_K.append(np.nanmedian(GEOS5FP_connection.Tmin_K(time_UTC=time_UTC, geometry=geometry)))
-            except Exception as e:
-                Tmin_K.append(np.nan)
-                logger.exception(e)
-
-        if not "SZA" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("calculating solar zenith angle")
-                SZA.append(calculate_SZA_from_datetime(time_UTC, lat, lon))
-            except Exception as e:
-                SZA.append(np.nan)
-                logger.exception(e)
-
-        if not "KG" in FLiES_inputs_from_calval_df.columns:
-            try:
-                logger.info("selecting Koppen Geiger climate classification")
-                KG.append(load_koppen_geiger(geometry=geometry)[1, 1][0][0])
-            except Exception as e:
-                KG.append(np.nan)
-                logger.exception(e)
-
-    if not "doy" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["doy"] = doy
-
-    if not "AOT" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["AOT"] = AOT
-    
-    if not "COT" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["COT"] = COT
-    
-    if not "vapor_gccm" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["vapor_gccm"] = vapor_gccm
-
-    if not "ozone_cm" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["ozone_cm"] = ozone_cm
-    
-    FLiES_inputs_df["SWin"] = SWin
-    FLiES_inputs_df["Tmin_K"] = Tmin_K 
+    Generates a DataFrame of FLiES inputs by retrieving atmospheric and static data.
     
-    if not "SZA" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["SZA"] = SZA
+    This is a simple wrapper around retrieve_FLiESANN_inputs that handles DataFrame
+    input/output and geometry parsing.
+
+    Parameters:
+    input_df (pd.DataFrame): A DataFrame containing the following columns:
+        - time_UTC (str or datetime): Time in UTC.
+        - geometry (str or shapely.geometry.Point) or (lat, lon): Spatial coordinates. 
+          If "geometry" is a string, it should be in WKT format (e.g., "POINT (lon lat)").
+        - albedo (float, optional): Surface albedo.
+        - COT (float, optional): Cloud optical thickness.
+        - AOT (float, optional): Aerosol optical thickness.
+        - vapor_gccm (float, optional): Water vapor in grams per cubic centimeter.
+        - ozone_cm (float, optional): Ozone concentration in centimeters.
+        - elevation_m (float, optional): Elevation in meters.
+        - SZA_deg (float, optional): Solar zenith angle in degrees.
+        - KG or KG_climate (str, optional): Köppen-Geiger climate classification.
+        - SWin_Wm2 (float, optional): Shortwave incoming solar radiation.
+        - day_of_year (float, optional): Day of year.
+        - hour_of_day (float, optional): Hour of day.
+    GEOS5FP_connection (GEOS5FP, optional): Connection object for GEOS-5 FP data.
+    NASADEM_connection (NASADEMConnection, optional): Connection object for NASADEM data.
+
+    Returns:
+    pd.DataFrame: A DataFrame with the same structure as the input, but with additional columns:
+        - albedo: Surface albedo array
+        - COT: Cloud optical thickness array
+        - AOT: Aerosol optical thickness array
+        - vapor_gccm: Water vapor array
+        - ozone_cm: Ozone concentration array
+        - elevation_m: Elevation in meters array
+        - elevation_km: Elevation in kilometers array
+        - KG_climate: Climate classification array
+        - SZA_deg: Solar zenith angle array
+        - SWin_Wm2: Shortwave incoming radiation array
+        - day_of_year: Day of year array
+        - atype: Aerosol type array
+        - ctype: Cloud type array
+
+    Raises:
+    KeyError: If required columns ("geometry" or "lat" and "lon") are missing.
+    """
+    def ensure_geometry(row):
+        if "geometry" in row:
+            if isinstance(row.geometry, str):
+                s = row.geometry.strip()
+                if s.startswith("POINT"):
+                    coords = s.replace("POINT", "").replace("(", "").replace(")", "").strip().split()
+                    return Point(float(coords[0]), float(coords[1]))
+                elif "," in s:
+                    coords = [float(c) for c in s.split(",")]
+                    return Point(coords[0], coords[1])
+                else:
+                    coords = [float(c) for c in s.split()]
+                    return Point(coords[0], coords[1])
+        return row.geometry
+
+    logger.info("started generating FLiES input table")
+
+    # Ensure geometry column is properly formatted
+    input_df = input_df.copy()
+    input_df["geometry"] = input_df.apply(ensure_geometry, axis=1)
+
+    # Prepare output DataFrame
+    output_df = input_df.copy()
 
-    if not "KG" in FLiES_inputs_df.columns:
-        FLiES_inputs_df["KG"] = KG
+    # Prepare geometries
+    if "geometry" in input_df.columns:
+        geometries = MultiPoint([(geom.x, geom.y) for geom in input_df.geometry], crs=WGS84)
+    elif "lat" in input_df.columns and "lon" in input_df.columns:
+        geometries = MultiPoint([(lon, lat) for lon, lat in zip(input_df.lon, input_df.lat)], crs=WGS84)
+    else:
+        raise KeyError("Input DataFrame must contain either 'geometry' or both 'lat' and 'lon' columns.")
 
-    if "Ta" in FLiES_inputs_df and "Ta_C" not in FLiES_inputs_df:
-        FLiES_inputs_df.rename({"Ta": "Ta_C"}, inplace=True)
+    # Convert time column to datetime
+    times_UTC = pd.to_datetime(input_df.time_UTC)
 
-    return FLiES_inputs_df
+    logger.info(f"generating inputs for {len(input_df)} rows")
+
+    # Helper function to get column values or None if column doesn't exist
+    def get_column_or_none(df, col_name, default_col_name=None):
+        if col_name in df.columns:
+            return df[col_name].values
+        elif default_col_name and default_col_name in df.columns:
+            return df[default_col_name].values
+        else:
+            return None
+
+    # Retrieve all inputs at once using vectorized retrieve_FLiESANN_inputs call
+    FLiES_inputs = retrieve_FLiESANN_inputs(
+        geometry=geometries,
+        time_UTC=times_UTC,
+        albedo=get_column_or_none(input_df, "albedo"),
+        COT=get_column_or_none(input_df, "COT"),
+        AOT=get_column_or_none(input_df, "AOT"),
+        vapor_gccm=get_column_or_none(input_df, "vapor_gccm"),
+        ozone_cm=get_column_or_none(input_df, "ozone_cm"),
+        elevation_m=get_column_or_none(input_df, "elevation_m"),
+        SZA_deg=get_column_or_none(input_df, "SZA_deg", "SZA"),
+        KG_climate=get_column_or_none(input_df, "KG_climate", "KG"),
+        SWin_Wm2=get_column_or_none(input_df, "SWin_Wm2"),
+        day_of_year=get_column_or_none(input_df, "day_of_year"),
+        hour_of_day=get_column_or_none(input_df, "hour_of_day"),
+        GEOS5FP_connection=GEOS5FP_connection,
+        NASADEM_connection=NASADEM_connection
+    )
+
+    # Add retrieved inputs to the output DataFrame
+    for key, values in FLiES_inputs.items():
+        # Skip values with mismatched lengths
+        if hasattr(values, '__len__') and not isinstance(values, str):
+            if len(values) != len(output_df):
+                continue
+        output_df[key] = values
+
+    logger.info("completed generating FLiES input table")
+
+    return output_df