1717from rasterio .warp import calculate_default_transform , reproject , Resampling , transform_geom
1818from rasterio .features import bounds as geom_bounds
1919import warnings
20+ import logging
2021
2122from .utlis import *
2223from .interactS3 import *
2728from ..streamflowdata .forecasteddata import getNWMForecasteddata
2829from ..runFIM import runOWPHANDFIM
2930
31+ logging .getLogger ("rasterio" ).setLevel (logging .ERROR )
32+ logging .getLogger ("rasterio._env" ).setLevel (logging .ERROR )
33+ warnings .filterwarnings ("ignore" , category = UserWarning )
34+
3035
3136#GET LOW FIDELITY USING FIMSERVE
3237def get_LFFIM (huc_id , event_date = None , data = 'forecast' , forecast_range = None , forecast_date = None , sort_by = None ):
@@ -347,7 +352,7 @@ def clip_raster_inplace_to_boundary(
347352 persists outside the boundary (no "black dots"/artifacts outside the clip).
348353 """
349354 raster_path = Path (raster_path )
350- clipped_path = raster_path .with_name (f"{ raster_path .stem } _Clipped { raster_path .suffix } " )
355+ clipped_path = raster_path .with_name (f"{ raster_path .stem } _clipped { raster_path .suffix } " )
351356
352357 geoms , b_crs = _ensure_list_of_geoms_and_crs (boundary_geometry , boundary_crs = boundary_crs )
353358 if not geoms :
@@ -377,9 +382,6 @@ def clip_raster_inplace_to_boundary(
377382 all_touched = False
378383 )
379384
380- # Extra safety: ensure masked-out pixels are exactly nodata_value across all bands
381- if src .count == 1 :
382- out_image = np .where (out_image == out_image , out_image , nodata_value ) # no-op for safety
383385 # mask() already filled; the key is to force dtype + nodata consistency
384386 if out_image .dtype != src .dtypes [0 ]:
385387 out_image = out_image .astype (src .dtypes [0 ], copy = False )
@@ -436,7 +438,8 @@ def clip_all_forcings_if_boundary_overlaps(
436438 empty dict if clipping was skipped due to non-overlap.
437439 """
438440 forcing_dir = Path (forcing_dir )
439- tifs = sorted (forcing_dir .glob ("*.tif" ))
441+ all_tifs = sorted (forcing_dir .glob ("*.tif" ))
442+ tifs = [p for p in all_tifs if not p .stem .endswith ("_clipped" )]
440443
441444 if not tifs :
442445 raise FileNotFoundError (f"No .tif forcing rasters found in: { forcing_dir } " )
@@ -454,12 +457,12 @@ def clip_all_forcings_if_boundary_overlaps(
454457 non_overlapping .append (tif )
455458 continue
456459
457- # Transform boundary geometries into this forcing CRS before overlap check
458460 if b_crs is None :
459- b_crs_local = src . crs
461+ b_crs_local = CRS . from_user_input ( boundary_crs )
460462 else :
461463 b_crs_local = b_crs
462464
465+ # Transform boundary geometries into this forcing CRS before overlap check
463466 if src .crs != b_crs_local :
464467 geoms_in_src = [transform_geom (b_crs_local , src .crs , g , precision = 6 ) for g in geoms ]
465468 else :
@@ -495,6 +498,66 @@ def clip_all_forcings_if_boundary_overlaps(
495498
496499 return mapping
497500
501+
502+ def clip_fim_to_boundary (
503+ fim_raster_path : Path ,
504+ boundary_geometry ,
505+ boundary_crs : Union [str , dict ] = "EPSG:4326" ,
506+ nodata_value : int = 0
507+ ) -> Path :
508+ fim_raster_path = Path (fim_raster_path )
509+ clipped_path = fim_raster_path .with_name (f"{ fim_raster_path .stem } _clipped{ fim_raster_path .suffix } " )
510+
511+ geoms , b_crs = _ensure_list_of_geoms_and_crs (boundary_geometry , boundary_crs = boundary_crs )
512+ if not geoms :
513+ raise ValueError ("boundary_geometry is empty; cannot clip." )
514+
515+ with rasterio .open (fim_raster_path ) as src :
516+ src_crs = src .crs
517+ if src_crs is None :
518+ raise ValueError (f"Raster has no CRS: { fim_raster_path } " )
519+
520+ if b_crs is None :
521+ b_crs = src_crs
522+
523+ if src_crs != b_crs :
524+ geoms_in_src = [transform_geom (b_crs , src_crs , g , precision = 6 ) for g in geoms ]
525+ else :
526+ geoms_in_src = geoms
527+
528+ out_image , out_transform = mask (
529+ src ,
530+ geoms_in_src ,
531+ crop = True ,
532+ filled = True ,
533+ nodata = nodata_value ,
534+ all_touched = False
535+ )
536+
537+ if out_image .dtype != src .dtypes [0 ]:
538+ out_image = out_image .astype (src .dtypes [0 ], copy = False )
539+
540+ out_meta = src .meta .copy ()
541+ out_meta .update ({
542+ "driver" : "GTiff" ,
543+ "height" : out_image .shape [1 ],
544+ "width" : out_image .shape [2 ],
545+ "transform" : out_transform ,
546+ "crs" : src_crs ,
547+ "nodata" : nodata_value ,
548+ })
549+
550+ with rasterio .open (clipped_path , "w" , ** out_meta ) as dst :
551+ dst .write (out_image )
552+
553+ try :
554+ compress_tif_lzw (clipped_path )
555+ except Exception :
556+ pass
557+
558+ return clipped_path
559+
560+
498561# PREPROCESS THE OWP HAND BASED FIM FOR SM
499562def prepare_FORCINGs (
500563 huc_id ,
@@ -514,8 +577,9 @@ def prepare_FORCINGs(
514577
515578 # If here, some boundary is passed, If that boundary overlaps with all the forcings,
516579 forcing_dir = Path (f'./HUC{ huc_id } _forcings' )
517- forcing_suffix = ""
518580
581+ mapping = {}
582+ did_clip_forcings = False
519583 if clip_boundary is not None :
520584 print ("Boundary provided. Checking overlap with all forcings...\n " )
521585 mapping = clip_all_forcings_if_boundary_overlaps (
@@ -524,7 +588,7 @@ def prepare_FORCINGs(
524588 boundary_crs = clip_boundary_crs
525589 )
526590 if mapping :
527- forcing_suffix = "_clipped"
591+ did_clip_forcings = True
528592 print ("All forcing rasters clipped successfully.\n " )
529593 else :
530594 print ("Skipping forcing clipping due to non-overlap.\n " )
@@ -543,6 +607,9 @@ def prepare_FORCINGs(
543607 # Get the HUC8 boundary
544608 HUC_boundary = getHUC8BoundaryByID (huc_id )
545609
610+ lulc_original = forcing_dir / f"LULC_HUC{ huc_id } .tif"
611+ reference_dir = mapping .get (lulc_original , lulc_original )
612+
546613 for FIM in fim_files :
547614 out_dir = fim_dir / 'processing'
548615 out_dir .mkdir (parents = True , exist_ok = True )
@@ -562,10 +629,24 @@ def prepare_FORCINGs(
562629 mask_with_PWB (out_dir_binary , final_raster )
563630 compress_tif_lzw (final_raster )
564631
632+ # If boundary clipping happened for forcings, clip the FIM as well
633+ final_for_alignment = final_raster
634+ if did_clip_forcings and clip_boundary is not None :
635+ final_for_alignment = clip_fim_to_boundary (
636+ fim_raster_path = final_raster ,
637+ boundary_geometry = clip_boundary ,
638+ boundary_crs = clip_boundary_crs ,
639+ nodata_value = 0
640+ )
641+
565642 # Align final FIM raster with reference raster
566- reference_dir = forcing_dir / f'LULC_HUC{ huc_id } { forcing_suffix } .tif'
567- FIM_finaldir = forcing_dir / f'hand_{ FIM .stem } .tif'
568- align_raster (final_raster , reference_dir , FIM_finaldir )
643+ fim_name = FIM .stem
644+ if did_clip_forcings and clip_boundary is not None :
645+ FIM_finaldir = forcing_dir / f'hand_{ fim_name } _clipped.tif'
646+ else :
647+ FIM_finaldir = forcing_dir / f'hand_{ fim_name } .tif'
648+
649+ align_raster (final_for_alignment , reference_dir , FIM_finaldir )
569650
570651 # Clean up temporary FIM directory
571652 if cwd .exists () and cwd .is_dir ():
0 commit comments