streamline the surrogate modeling within fimserve; published new PyPi version

Author: supathdhitalGEO

build/lib/fimserve/FIMsubset/__init__.py

@@ -0,0 +1,7 @@
+import warnings
+
+warnings.simplefilter("ignore")
+
+from .xycoord import subsetFIM
+
+__all__ = ["subsetFIM"]

build/lib/fimserve/FIMsubset/shpsubset.py

@@ -0,0 +1,61 @@
+import os
+import glob
+import rasterio
+from rasterio.mask import mask
+from pathlib import Path
+import geopandas as gpd
+
+from ..datadownload import setup_directories
+
+
+def checkSHP(input_file):
+    if isinstance(input_file, Path):
+        input_file = str(input_file)
+
+    if input_file.endswith(".shp"):
+        gdf = gpd.read_file(input_file)
+    elif input_file.endswith(".gpkg"):
+        gdf = gpd.read_file(input_file)
+    elif input_file.endswith(".geojson"):
+        gdf = gpd.read_file(input_file)
+    elif input_file.endswith(".kml"):
+        gdf = gpd.read_file(input_file)
+    else:
+        raise ValueError(
+            "Unsupported file format. Please provide a shapefile GeoPackage, GeoJSON, or KML."
+        )
+    return gdf
+
+
+def clipFIMforboundary(inundation_raster, shapefile_geom, output_directory, huc):
+    with rasterio.open(inundation_raster) as src:
+        geometries = shapefile_geom.geometry.tolist()
+        if src.crs != shapefile_geom.crs:
+            print("CRS mismatch. Reprojecting geometry to match raster.")
+            geometries = [geom.to_crs(src.crs) for geom in geometries]
+
+        out_image, out_transform = mask(src, geometries, crop=True)
+        out_meta = src.meta.copy()
+        out_meta.update(
+            {
+                "driver": "GTiff",
+                "count": 1,
+                "crs": src.crs,
+                "transform": out_transform,
+                "width": out_image.shape[2],
+                "height": out_image.shape[1],
+            }
+        )
+
+        output_directory = os.path.join(output_directory, "subsetFIM")
+        if not os.path.exists(output_directory):
+            os.makedirs(output_directory)
+
+        file_basename = os.path.basename(inundation_raster).split("_")[0]
+        output_file = os.path.join(output_directory, f"{file_basename}_subsetFIM.tif")
+
+        if not os.path.exists(output_directory):
+            os.makedirs(output_directory)
+        with rasterio.open(output_file, "w", **out_meta) as dest:
+            dest.write(out_image)
+        print(f"Clipped raster saved to {output_file}")

build/lib/fimserve/FIMsubset/xycoord.py

@@ -0,0 +1,102 @@
+import os
+import glob
+import geopandas as gpd
+from shapely.geometry import Point
+from pyproj import CRS
+import rasterio
+from rasterio.mask import mask
+
+from ..datadownload import setup_directories
+from .shpsubset import checkSHP, clipFIMforboundary
+
+
+def checkifWGS(x, y):
+    """Check if coordinates are in WGS84 (EPSG:4326) range."""
+    if -180 <= x <= 180 and -90 <= y <= 90:
+        return True
+    else:
+        return False
+
+
+def reproject_coordinates(x, y, target_crs="EPSG:5070"):
+    """Reproject coordinates to target CRS."""
+    if checkifWGS(x, y):
+        from pyproj import CRS, Transformer
+
+        crs_wgs84 = CRS("EPSG:4326")
+        crs_target = CRS(target_crs)
+        transformer = Transformer.from_crs(crs_wgs84, crs_target, always_xy=True)
+        x_new, y_new = transformer.transform(x, y)
+        return x_new, y_new
+    else:
+        return x, y
+
+
+def clipFIM(inundation_raster, shapefile_geom, output_directory, huc):
+    with rasterio.open(inundation_raster) as src:
+        out_image, out_transform = mask(src, [shapefile_geom], crop=True)
+        out_meta = src.meta.copy()
+        out_meta.update(
+            {
+                "driver": "GTiff",
+                "count": 1,
+                "crs": src.crs,
+                "transform": out_transform,
+                "width": out_image.shape[2],
+                "height": out_image.shape[1],
+            }
+        )
+        file_basename = os.path.basename(inundation_raster).split("_")[0]
+        output_file = os.path.join(output_directory, f"{file_basename}_subsetFIM.tif")
+
+        if not os.path.exists(output_directory):
+            os.makedirs(output_directory)
+
+        with rasterio.open(output_file, "w", **out_meta) as dest:
+            dest.write(out_image)
+
+        print(f"Clipped raster saved to {output_file}")
+        return out_image, out_transform
+
+
+def withininWatershed(x, y, geopackage_path, inundation_raster, huc, output_directory):
+    point = Point(x, y)
+    watersheds = gpd.read_file(geopackage_path)
+
+    # Watershed containing the point
+    watershed_containing_point = None
+    for _, watershed in watersheds.iterrows():
+        if watershed["geometry"].contains(point):
+            watershed_containing_point = watershed
+            break
+    output_directory = os.path.join(output_directory, "subsetFIM")
+    if not os.path.exists(output_directory):
+        os.makedirs(output_directory)
+
+    if watershed_containing_point is None:
+        raise ValueError(f"Point is not within any watershed boundary for {huc}")
+    return clipFIM(
+        inundation_raster, watershed_containing_point["geometry"], output_directory, huc
+    )
+
+
+def subsetFIM(location, huc, method):
+    code_dir, data_dir, output_dir = setup_directories()
+    gpkg_path = os.path.join(
+        output_dir, f"flood_{huc}", huc, "nwm_catchments_proj_subset.gpkg"
+    )
+    out_dir = os.path.join(output_dir, f"flood_{huc}", f"{huc}_inundation")
+    inundation_rasters = os.path.join(out_dir, "*_inundation.tif")
+    print(inundation_rasters)
+    files = glob.glob(inundation_rasters)
+    if method == "xy":
+        x, y = location
+        print("point", x, y)
+        x, y = reproject_coordinates(x, y)
+        for file in files:
+            print(f"Clipping {file} to watershed containing point {x}, {y}")
+            withininWatershed(x, y, gpkg_path, file, huc, out_dir)
+    elif method == "boundary":
+        shapefile = checkSHP(location)
+        for file in files:
+            clipFIMforboundary(file, shapefile, out_dir, huc)

build/lib/fimserve/__init__.py

@@ -0,0 +1,84 @@
+import warnings
+
+from fimserve.fimevaluation.fims_setup import fim_lookup
+
+warnings.simplefilter("ignore")
+
+from .datadownload import DownloadHUC8
+from .streamflowdata.nwmretrospective import getNWMretrospectivedata
+from .runFIM import runOWPHANDFIM
+
+from .streamflowdata.forecasteddata import getNWMForecasteddata
+from .streamflowdata.geoglows import getGEOGLOWSstreamflow
+
+# plots
+from .plot.nwmfid import plotNWMStreamflow
+from .streamflowdata.usgsdata import getUSGSsitedata
+from .plot.comparestreamflow import CompareNWMnUSGSStreamflow
+from .plot.usgs import plotUSGSStreamflow
+from .plot.src import plotSRC
+
+# For table
+from .plot.usgsandfid import GetUSGSIDandCorrFID
+
+# subsetting
+from .FIMsubset.xycoord import subsetFIM
+
+# Fim visualization
+from .vizualizationFIM import vizualizeFIM
+
+
+# Statistics
+from .statistics.calculatestatistics import CalculateStatistics
+
+# For intersected HUC8 boundary
+from .intersectedHUC import getIntersectedHUC8ID
+
+
+# evaluation of FIM
+from .fimevaluation.fims_setup import FIMService, fim_lookup
+from .fimevaluation.run_fimeval import run_evaluation
+
+
+#Enhancement using surrogate models [Importing those only if they are called]
+def prepare_FORCINGs(*args, **kwargs):
+    from .enhancement_withSM.preprocessFIM import prepare_FORCINGs as _impl
+    return _impl(*args, **kwargs)
+
+def predict_SM(*args, **kwargs):
+    from .enhancement_withSM.SM_prediction import predict_SM as _impl
+    return _impl(*args, **kwargs)
+
+def get_building_exposure(*args, **kwargs):
+    from .enhancement_withSM.building_exposure import get_building_exposure as _impl
+    return _impl(*args, **kwargs)
+
+def get_population_exposure(*args, **kwargs):
+    from .enhancement_withSM.pop_exposure import get_population_exposure as _impl
+    return _impl(*args, **kwargs)
+
+
+__all__ = [
+    "DownloadHUC8",
+    "getNWMRetrospectivedata",
+    "runOWPHANDFIM",
+    "getNWMForecasteddata",
+    "getGEOGLOWSstreamflow",
+    "plotNWMStreamflow",
+    "getUSGSsitedata",
+    "CompareNWMnUSGSStreamflow",
+    "plotUSGSStreamflow",
+    "plotSRC",
+    "GetUSGSIDandCorrFID",
+    "subsetFIM",
+    "vizualizeFIM",
+    "CalculateStatistics",
+    "getIntersectedHUC8ID",
+    "FIMService",
+    "fim_lookup",
+    "run_evaluation",
+    "prepare_FORCINGs",
+    "predict_SM",
+    "get_building_exposure",
+    "get_population_exposure",
+]