Code changes to EPSG:5070

data/nld/preprocess_levee_protected_areas.py

@@ -9,13 +9,14 @@
 from datetime import datetime
 
 sys.path.append('/foss_fim/src')
+from utils.shared_variables import DEFAULT_FIM_PROJECTION_CRS
 import utils.shared_functions as sf
 from utils.shared_functions import FIM_Helpers as fh
 
 def preprocess_levee_protected_areas(source_file_name_and_path,
                                      target_output_folder_path = '',
                                      target_output_filename = 'levee_protected_areas.gpkg',
-                                     projection_crs_id = 'EPSG:5070',
+                                     projection_crs_id = DEFAULT_FIM_PROJECTION_CRS,
                                      overwrite = False):
 
     '''
@@ -77,7 +78,7 @@ def preprocess_levee_protected_areas(source_file_name_and_path,
 
     # projection crs
     if (projection_crs_id is None) or (projection_crs_id == ''):
-        projection_crs_id = 'EPSG:5070'
+        projection_crs_id = DEFAULT_FIM_PROJECTION_CRS
 
     # -------------------
     # setup logs

data/usgs/acquire_and_preprocess_3dep_dems.py

@@ -143,7 +143,7 @@ def __download_usgs_dems(extent_files, output_folder_path, number_of_jobs, retry
     base_cmd += ' -cutline {2} -crop_to_cutline -ot Float32 -r bilinear'
     base_cmd += ' -of "GTiff" -overwrite -co "BLOCKXSIZE=256" -co "BLOCKYSIZE=256"'
     base_cmd += ' -co "TILED=YES" -co "COMPRESS=LZW" -co "BIGTIFF=YES" -tr 10 10'
-    base_cmd += ' -t_srs EPSG:5070 -cblend 6'
+    base_cmd += ' -t_srs {3} -cblend 6'
 
     with ProcessPoolExecutor(max_workers=number_of_jobs) as executor:
 
@@ -212,7 +212,7 @@ def download_usgs_dem_file(extent_file,
                 base_cmd += ' -cutline {2} -crop_to_cutline -ot Float32 -r bilinear'
                 base_cmd +=  ' -of "GTiff" -overwrite -co "BLOCKXSIZE=256" -co "BLOCKYSIZE=256"'
                 base_cmd +=  ' -co "TILED=YES" -co "COMPRESS=LZW" -co "BIGTIFF=YES" -tr 10 10'
-                base_cmd +=  ' -t_srs EPSG:5070 -cblend 6'
+                base_cmd +=  ' -t_srs {3} -cblend 6'
         - retry (bool)
              If True, and the file exists (and is over 0k), downloading will be skipped.
         
@@ -246,7 +246,8 @@ def download_usgs_dem_file(extent_file,
 
     cmd = base_cmd.format(download_url,
                             target_path_raw,
-                            extent_file)
+                            extent_file,
+                            sv.DEFAULT_FIM_PROJECTION_CRS)
     #PREP_PROJECTION_EPSG
     #fh.vprint(f"cmd is {cmd}", self.is_verbose, True)
     #print(f"cmd is {cmd}")

gms_run_post_processing.sh

@@ -141,9 +141,9 @@ if [ "$src_adjust_spatial" = "True" ]; then
         export CALIBRATION_DB_PASS=$CALIBRATION_DB_PASS
         echo "Populate PostgrSQL database with benchmark FIM extent points and HUC attributes (the calibration database)"
         echo "Loading HUC Data"
-        time ogr2ogr -overwrite -nln hucs -a_srs ESRI:102039 -f PostgreSQL PG:"host=$CALIBRATION_DB_HOST dbname=$CALIBRATION_DB_NAME user=$CALIBRATION_DB_USER_NAME password=$CALIBRATION_DB_PASS" $inputDataDir/wbd/WBD_National.gpkg WBDHU8
+        time ogr2ogr -overwrite -nln hucs -t_srs EPSG:5070 -f PostgreSQL PG:"host=$CALIBRATION_DB_HOST dbname=$CALIBRATION_DB_NAME user=$CALIBRATION_DB_USER_NAME password=$CALIBRATION_DB_PASS" $inputDataDir/wbd/WBD_National.gpkg WBDHU8
         echo "Loading Point Data"
-        time ogr2ogr -overwrite -f PostgreSQL PG:"host=$CALIBRATION_DB_HOST dbname=$CALIBRATION_DB_NAME user=$CALIBRATION_DB_USER_NAME password=$CALIBRATION_DB_PASS" $fim_obs_pnt_data usgs_nws_benchmark_points -nln points
+        time ogr2ogr -overwrite -t_srs EPSG:5070 -f PostgreSQL PG:"host=$CALIBRATION_DB_HOST dbname=$CALIBRATION_DB_NAME user=$CALIBRATION_DB_USER_NAME password=$CALIBRATION_DB_PASS" $fim_obs_pnt_data usgs_nws_benchmark_points -nln points
     fi
 fi
 

src/clip_vectors_to_wbd.py

@@ -6,6 +6,7 @@
 import rasterio as rio
 
 from shapely.geometry import MultiPolygon,Polygon
+from utils.shared_variables import DEFAULT_FIM_PROJECTION_CRS
 from utils.shared_functions import getDriver, mem_profile
 
 @mem_profile
@@ -83,7 +84,7 @@ def subset_vector_layers(subset_nwm_lakes,
     print("Subsetting Levee Protected Areas", flush=True)
     levee_protected_areas = gpd.read_file(levee_protected_areas, mask=wbd_buffer)
     if not levee_protected_areas.empty:
-        levee_protected_areas = levee_protected_areas.to_crs('ESRI:102039')
+        # levee_protected_areas = levee_protected_areas.to_crs(DEFAULT_FIM_PROJECTION_CRS)
         levee_protected_areas.to_file(subset_levee_protected_areas, driver = getDriver
                                       (subset_levee_protected_areas), index=False)
     del levee_protected_areas

src/src_adjust_spatial_obs.py

@@ -18,7 +18,7 @@
 from geopandas.tools import sjoin
 from multiprocessing import Pool
 from src_roughness_optimization import update_rating_curve
-from utils.shared_variables import DOWNSTREAM_THRESHOLD, ROUGHNESS_MIN_THRESH, ROUGHNESS_MAX_THRESH
+from utils.shared_variables import DOWNSTREAM_THRESHOLD, ROUGHNESS_MIN_THRESH, ROUGHNESS_MAX_THRESH, DEFAULT_FIM_PROJECTION_CRS
 
 #import variables from .env file
 load_dotenv()
@@ -52,7 +52,7 @@
 def process_points(args):
 
     '''
-    The funciton ingests geodataframe and attributes the point data with its hydroid and HAND values before passing a dataframe to the src_roughness_optimization.py workflow
+    The function ingests geodataframe and attributes the point data with its hydroid and HAND values before passing a dataframe to the src_roughness_optimization.py workflow
 
     Processing
     - Extract x,y coordinates from geometry
@@ -74,21 +74,12 @@ def process_points(args):
     ## Define coords variable to be used in point raster value attribution.
     coords = [(x,y) for x, y in zip(water_edge_df.X, water_edge_df.Y)]
 
+    water_edge_df.to_crs(DEFAULT_FIM_PROJECTION_CRS)
+
     ## Use point geometry to determine HAND raster pixel values.
-    hand_src = rasterio.open(hand_path)
-    hand_crs = hand_src.crs
-    water_edge_df.to_crs(hand_crs)  # Reproject geodataframe to match hand_src. Should be the same, but this is a double check.
-    water_edge_df['hand'] = [h[0] for h in hand_src.sample(coords)]
-    hand_src.close()
-    del hand_src, hand_crs,
-
-    ## Use point geometry to determine catchment raster pixel values.
-    catchments_src = rasterio.open(catchments_path)
-    catchments_crs = catchments_src.crs
-    water_edge_df.to_crs(catchments_crs)
-    water_edge_df['hydroid'] = [c[0] for c in catchments_src.sample(coords)]
-    catchments_src.close()
-    del catchments_src, catchments_crs
+    with rasterio.open(hand_path) as hand_src, rasterio.open(catchments_path) as catchments_src:
+        water_edge_df['hand'] = [h[0] for h in hand_src.sample(coords)]
+        water_edge_df['hydroid'] = [c[0] for c in catchments_src.sample(coords)]
 
     water_edge_df = water_edge_df[(water_edge_df['hydroid'].notnull()) & (water_edge_df['hand'] > 0) & (water_edge_df['hydroid'] > 0)]
 
@@ -156,10 +147,10 @@ def find_points_in_huc(huc_id, conn):
     JOIN hucs H ON ST_Contains(H.geom, P.geom)
     WHERE H.huc8 = %s """
 
-    # Need to hard code the CRS to use EPSG:5070 instead of the default ESRI:102039 (gdal pyproj throws an error with crs 102039)
+    # Use EPSG:5070 instead of the default ESRI:102039 (gdal pyproj throws an error with crs 102039)
     # Appears that EPSG:5070 is functionally equivalent to ESRI:102039: https://gis.stackexchange.com/questions/329123/crs-interpretation-in-qgis
     water_edge_df = gpd.GeoDataFrame.from_postgis(huc_pt_query, con=conn, 
-                                                  params=[huc_id], crs="EPSG:5070",
+                                                  params=[huc_id], crs=DEFAULT_FIM_PROJECTION_CRS,
                                                   parse_dates=['coll_time'])
     water_edge_df = water_edge_df.drop(columns=['st_x','st_y'])
 

src/utils/shared_variables.py

@@ -5,7 +5,7 @@
 
 # -- Projections: EPSG 5070-- #
 PREP_PROJECTION = 'PROJCRS["NAD83 / Conus Albers",BASEGEOGCRS["NAD83",DATUM["North American Datum 1983",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4269]],CONVERSION["Conus Albers",METHOD["Albers Equal Area",ID["EPSG",9822]],PARAMETER["Latitude of false origin",23,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",-96,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",29.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",45.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Data analysis and small scale data presentation for contiguous lower 48 states."],AREA["United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming."],BBOX[24.41,-124.79,49.38,-66.91]],ID["EPSG",5070]]' 
-DEFAULT_FIM_PROJECTION_CRS = 'ESPG:5070'  
+DEFAULT_FIM_PROJECTION_CRS = 'EPSG:5070'  
 PREP_CRS = CRS(PREP_PROJECTION)
 VIZ_PROJECTION ='PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]'
 

tools/inundate_nation.py

@@ -16,7 +16,7 @@
 from gms_tools.mosaic_inundation import Mosaic_inundation
 from gms_tools.inundate_gms import Inundate_gms
 from inundation import inundate
-from tools_shared_variables import elev_raster_ndv
+from utils.shared_variables import elev_raster_ndv, PREP_PROJECTION
 from utils.shared_functions import FIM_Helpers as fh
 
 
@@ -25,7 +25,7 @@
 #INPUTS_DIR = r'/data/inputs'
 #OUTPUT_BOOL_PARENT_DIR = '/data/inundation_review/inundate_nation/bool_temp/'
 #DEFAULT_OUTPUT_DIR = '/data/inundation_review/inundate_nation/mosaic_output/'
-PREP_PROJECTION = 'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010042,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'
+# PREP_PROJECTION = 'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010042,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'
 
 
 def inundate_nation(fim_run_dir, output_dir, magnitude_key, flow_file, inc_mosaic, job_number):