Anthro updates

lib/commons/rscommons/moving_window.py

@@ -47,14 +47,14 @@ def moving_window_dgo_ids(igo: str, dgo: str, level_paths: list, distance: dict)
 
     with GeopackageLayer(igo) as lyr_igo, GeopackageLayer(dgo) as lyr_dgo:
         for level_path in level_paths:
-            for feat_seg_pt, *_, in lyr_igo.iterate_features(f'Finding windows on {level_path}', attribute_filter=f"LevelPathI = {level_path}"):
+            for feat_seg_pt, *_, in lyr_igo.iterate_features(f'Finding windows on {level_path}', attribute_filter=f"level_path = {level_path}"):
                 window_distance = distance[str(feat_seg_pt.GetField('stream_size'))]
                 dist = feat_seg_pt.GetField('seg_distance')
                 min_dist = dist - 0.5 * window_distance
                 max_dist = dist + 0.5 * window_distance
 
                 dgoids = []
-                for feat_seg_poly, *_ in lyr_dgo.iterate_features(attribute_filter=f"LevelPathI = {int(level_path)} and seg_distance >= {int(min_dist)} and seg_distance <= {int(max_dist)}"):
+                for feat_seg_poly, *_ in lyr_dgo.iterate_features(attribute_filter=f"level_path = {int(level_path)} and seg_distance >= {int(min_dist)} and seg_distance <= {int(max_dist)}"):
                     dgoids.append(feat_seg_poly.GetFID())
                 windows[feat_seg_pt.GetFID()] = dgoids
 

packages/anthro/.vscode/launch.json

@@ -26,7 +26,7 @@
                 "{env:DATA_ROOT}/rs_context/${input:HUC}/topography/dem_hillshade.tif",
                 "{env:DATA_ROOT}/vbet/${input:HUC}/outputs/vbet.gpkg/vbet_igos",
                 "{env:DATA_ROOT}/vbet/${input:HUC}/intermediates/vbet_intermediates.gpkg/vbet_dgos",
-                "{env:DATA_ROOT}/rs_context/${input:HUC}/hydrology/hydro_derivatives.gpkg/network_intersected_300m",
+                "{env:DATA_ROOT}/rs_context/${input:HUC}/hydrology/hydro_derivatives.gpkg/network_intersected",
                 "{env:DATA_ROOT}/vbet/${input:HUC}/outputs/vbet.gpkg/vbet_full",
                 "{env:DATA_ROOT}/rs_context/${input:HUC}/ownership/ownership.shp",
                 "{env:DATA_ROOT}/rs_context/${input:HUC}/transportation/canals.shp",

packages/anthro/anthro/__version__.py

@@ -1 +1 @@
-__version__ = "0.2.0"
+__version__ = "0.2.1"

packages/anthro/anthro/anthro.py

@@ -134,7 +134,8 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
     # Create the output feature class fields. Only those listed here will get copied from the source.
     with GeopackageLayer(outputs_gpkg_path, layer_name=LayerTypes['OUTPUTS'].sub_layers['ANTHRO_GEOM_POINTS'].rel_path, delete_dataset=True) as out_lyr:
         out_lyr.create_layer(ogr.wkbMultiPoint, epsg=cfg.OUTPUT_EPSG, options=['FID=IGOID'], fields={
-            'LevelPathI': ogr.OFTReal,
+            'FCode': ogr.OFTInteger,
+            'level_path': ogr.OFTReal,
             'seg_distance': ogr.OFTReal,
             'stream_size': ogr.OFTInteger,
             'window_size': ogr.OFTReal,
@@ -144,7 +145,8 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
 
     with GeopackageLayer(outputs_gpkg_path, layer_name=LayerTypes['OUTPUTS'].sub_layers['ANTHRO_GEOM_DGOS'].rel_path, write=True) as out_lyr:
         out_lyr.create_layer(ogr.wkbMultiPolygon, epsg=cfg.OUTPUT_EPSG, options=['FID=DGOID'], fields={
-            'LevelPathI': ogr.OFTReal,
+            'FCode': ogr.OFTInteger,
+            'level_path': ogr.OFTReal,
             'seg_distance': ogr.OFTReal,
             'centerline_length': ogr.OFTReal,
             'segment_area': ogr.OFTReal
@@ -176,8 +178,8 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
 
     with SQLiteCon(outputs_gpkg_path) as database:
         database.curs.execute('INSERT INTO ReachAttributes (ReachID, Orig_DA, iGeo_DA, ReachCode, WatershedID, StreamName, NHDPlusID) SELECT ReachID, TotDASqKm, DivDASqKm, FCode, WatershedID, GNIS_NAME, NHDPlusID FROM ReachGeometry')
-        database.curs.execute('INSERT INTO IGOAttributes (IGOID, LevelPathI, seg_distance, stream_size) SELECT IGOID, LevelPathI, seg_distance, stream_size FROM IGOGeometry')
-        database.curs.execute('INSERT INTO DGOAttributes (DGOID, LevelPathI, seg_distance, segment_area, centerline_length) SELECT DGOID, LevelPathI, seg_distance, segment_area, centerline_length FROM DGOGeometry')
+        database.curs.execute('INSERT INTO IGOAttributes (IGOID, FCode, level_path, seg_distance, stream_size) SELECT IGOID, FCode, level_path, seg_distance, stream_size FROM IGOGeometry')
+        database.curs.execute('INSERT INTO DGOAttributes (DGOID, FCode, level_path, seg_distance, segment_area, centerline_length) SELECT DGOID, FCode, level_path, seg_distance, segment_area, centerline_length FROM DGOGeometry')
 
         # Register vwReaches as a feature layer as well as its geometry column
         database.curs.execute("""INSERT INTO gpkg_contents (table_name, data_type, identifier, min_x, min_y, max_x, max_y, srs_id)
@@ -198,20 +200,20 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
         database.curs.execute("""INSERT INTO gpkg_geometry_columns (table_name, column_name, geometry_type_name, srs_id, z, m)
             SELECT 'vwDgos', column_name, geometry_type_name, srs_id, z, m FROM gpkg_geometry_columns WHERE table_name = 'DGOGeometry'""")
 
-        database.conn.execute('CREATE INDEX ix_igo_levelpath on IGOGeometry(LevelPathI)')
+        database.conn.execute('CREATE INDEX ix_igo_levelpath on IGOGeometry(level_path)')
         database.conn.execute('CREATE INDEX ix_igo_segdist on IGOGeometry(seg_distance)')
         database.conn.execute('CREATE INDEX ix_igo_size on IGOGeometry(stream_size)')
-        database.conn.execute('CREATE INDEX ix_dgo_levelpath on DGOGeometry(LevelPathI)')
+        database.conn.execute('CREATE INDEX ix_dgo_levelpath on DGOGeometry(level_path)')
         database.conn.execute('CREATE INDEX ix_dgo_segdist on DGOGeometry(seg_distance)')
 
         database.conn.commit()
 
-        database.curs.execute('SELECT DISTINCT LevelPathI FROM IGOGeometry')
+        database.curs.execute('SELECT DISTINCT level_path FROM IGOGeometry')
         levelps = database.curs.fetchall()
-        levelpathsin = [lp['LevelPathI'] for lp in levelps]
+        levelpathsin = [lp['level_path'] for lp in levelps]
 
     with SQLiteCon(inputs_gpkg_path) as db:
-        db.conn.execute('CREATE INDEX ix_dgo_levelpath on dgo(LevelPathI)')
+        db.conn.execute('CREATE INDEX ix_dgo_levelpath on dgo(level_path)')
         db.conn.execute('CREATE INDEX ix_dgo_segdist on dgo(seg_distance)')
         db.conn.commit()
 

packages/anthro/anthro/utils/igo_infrastructure.py

@@ -4,8 +4,6 @@
 """
 import os
 import sqlite3
-import json
-import numpy as np
 from osgeo import ogr
 from rscommons import Logger, get_shp_or_gpkg, GeopackageLayer
 from rscommons.classes.vector_base import VectorBase, get_utm_zone_epsg
@@ -32,15 +30,17 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
         if ftr is not None:
             long = ftr.GetGeometryRef().GetEnvelope()[0]
             epsg_proj = get_utm_zone_epsg(long)
+            ref_layer = road
         else:
             epsg_proj = None
 
     if not epsg_proj:
-        with get_shp_or_gpkg(road) as inref:
+        with get_shp_or_gpkg(rail) as inref:
             ftr = inref.ogr_layer.GetFeature(1)
             if ftr is not None:
                 long = ftr.GetGeometryRef().GetEnvelope()[0]
                 epsg_proj = get_utm_zone_epsg(long)
+                ref_layer = rail
             else:
                 epsg_proj = None
 
@@ -50,6 +50,7 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
             if ftr is not None:
                 long = ftr.GetGeometryRef().GetEnvelope()[0]
                 epsg_proj = get_utm_zone_epsg(long)
+                ref_layer = canal
             else:
                 epsg_proj = None
 
@@ -58,15 +59,16 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
 
         return
 
-    with get_shp_or_gpkg(road) as reflyr:
+    with get_shp_or_gpkg(ref_layer) as reflyr:
         sref, transform = reflyr.get_transform_from_epsg(reflyr.spatial_ref, epsg_proj)
 
     attribs = {}
     gpkg_driver = ogr.GetDriverByName('GPKG')
     with GeopackageLayer(out_gpkg_path, 'DGOGeometry') as dgo_lyr:
         for dgo_ftr, *_ in dgo_lyr.iterate_features('initializing Attributes'):
             fid = dgo_ftr.GetFID()
-            attribs[fid] = {'Road_len': 0, 'Rail_len': 0, 'Canal_len': 0, 'Roadx_ct': 0, 'DivPts_ct': 0}
+            attribs[fid] = {'Road_len': 0, 'Rail_len': 0, 'Canal_len': 0, 'RoadX_ct': 0, 'DivPts_ct': 0, 'Road_prim_len': 0, 'Road_sec_len': 0, 'Road_4wd_len': 0}
+        dgo_ftr = None
         for dataset, label in in_data.items():
             log.info(f'Calculating metrics for dataset: {label}')
             dsrc = gpkg_driver.Open(os.path.dirname(dataset))
@@ -108,15 +110,49 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
                     else:
                         attribs[dgoid][lb1] = 1
 
+            dgo_ftr = None
+
+    log.info('Calculating road metrics by road type')
+    with GeopackageLayer(out_gpkg_path, 'DGOGeometry') as dgo_lyr, GeopackageLayer(road) as road_lyr:
+        if road_lyr.ogr_layer.GetFeatureCount() == 0:
+            log.info('No road features, skipping road length calculation')
+        else:
+            for dgo_ftr, *_ in dgo_lyr.iterate_features():
+                dgoid = dgo_ftr.GetFID()
+                dgo_geom = dgo_ftr.GetGeometryRef()
+                for road_ftr, *_ in road_lyr.iterate_features(clip_shape=dgo_geom, attribute_filter='tnmfrc in (1, 2, 3, 5)'):
+                    road_geom = road_ftr.GetGeometryRef()
+                    if road_geom.Intersects(dgo_geom):
+                        road_clip = road_geom.Intersection(dgo_geom)
+                        road_shapely = VectorBase.ogr2shapely(road_clip, transform=transform)
+                        attribs[dgoid]['Road_prim_len'] += road_shapely.length
+                road_ftr = None
+                for road_ftr, *_ in road_lyr.iterate_features(clip_shape=dgo_geom, attribute_filter='tnmfrc = 4'):
+                    road_geom = road_ftr.GetGeometryRef()
+                    if road_geom.Intersects(dgo_geom):
+                        road_clip = road_geom.Intersection(dgo_geom)
+                        road_shapely = VectorBase.ogr2shapely(road_clip, transform=transform)
+                        attribs[dgoid]['Road_sec_len'] += road_shapely.length
+                road_ftr = None
+                for road_ftr, *_ in road_lyr.iterate_features(clip_shape=dgo_geom, attribute_filter='tnmfrc = 6'):
+                    road_geom = road_ftr.GetGeometryRef()
+                    if road_geom.Intersects(dgo_geom):
+                        road_clip = road_geom.Intersection(dgo_geom)
+                        road_shapely = VectorBase.ogr2shapely(road_clip, transform=transform)
+                        attribs[dgoid]['Road_4wd_len'] += road_shapely.length
+
     conn = sqlite3.connect(out_gpkg_path)
     curs = conn.cursor()
 
     for dgoid, vals in attribs.items():
         curs.execute(f'UPDATE DGOAttributes SET Road_len = {vals["Road_len"]} WHERE DGOID = {dgoid}')
         curs.execute(f'UPDATE DGOAttributes SET Rail_len = {vals["Rail_len"]} WHERE DGOID = {dgoid}')
         curs.execute(f'UPDATE DGOAttributes SET Canal_len = {vals["Canal_len"]} WHERE DGOID = {dgoid}')
-        curs.execute(f'UPDATE DGOAttributes SET RoadX_ct = {vals["Roadx_ct"]} WHERE DGOID = {dgoid}')
+        curs.execute(f'UPDATE DGOAttributes SET RoadX_ct = {vals["RoadX_ct"]} WHERE DGOID = {dgoid}')
         curs.execute(f'UPDATE DGOAttributes SET DivPts_ct = {vals["DivPts_ct"]} WHERE DGOID = {dgoid}')
+        curs.execute(f'UPDATE DGOAttributes SET Road_prim_len = {vals["Road_prim_len"]} WHERE DGOID = {dgoid}')
+        curs.execute(f'UPDATE DGOAttributes SET Road_sec_len = {vals["Road_sec_len"]} WHERE DGOID = {dgoid}')
+        curs.execute(f'UPDATE DGOAttributes SET Road_4wd_len = {vals["Road_4wd_len"]} WHERE DGOID = {dgoid}')
 
     # summarize metrics from DGOs to IGOs using moving windows
     for igoid, dgoids in windows.items():
@@ -125,39 +161,57 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
         canal_len = 0
         roadx_ct = 0
         divpts_ct = 0
-        window_area = 0
+        road_prim_len = 0
+        road_sec_len = 0
+        road_4wd_len = 0
+        window_len = 0
 
         for dgoid in dgoids:
             road_len += curs.execute(f'SELECT Road_len FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
             rail_len += curs.execute(f'SELECT Rail_len FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
             canal_len += curs.execute(f'SELECT Canal_len FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
             roadx_ct += curs.execute(f'SELECT RoadX_ct FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
             divpts_ct += curs.execute(f'SELECT DivPts_ct FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
-            window_area += curs.execute(f'SELECT segment_area FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
+            road_prim_len += curs.execute(f'SELECT Road_prim_len FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
+            road_sec_len += curs.execute(f'SELECT Road_sec_len FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
+            road_4wd_len += curs.execute(f'SELECT Road_4wd_len FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
+            window_len += curs.execute(f'SELECT centerline_length FROM DGOAttributes WHERE DGOID = {dgoid}').fetchone()[0]
 
-        if window_area == 0:
+        if window_len == 0:
             road_dens = 0
             rail_dens = 0
             canal_dens = 0
             roadx_dens = 0
             divpts_dens = 0
+            road_prim_dens = 0
+            road_sec_dens = 0
+            road_4wd_dens = 0
         else:
-            road_dens = road_len / window_area
-            rail_dens = rail_len / window_area
-            canal_dens = canal_len / window_area
-            roadx_dens = roadx_ct / window_area
-            divpts_dens = divpts_ct / window_area
+            road_dens = road_len / window_len
+            rail_dens = rail_len / window_len
+            canal_dens = canal_len / window_len
+            roadx_dens = roadx_ct / window_len
+            divpts_dens = divpts_ct / window_len
+            road_prim_dens = road_prim_len / window_len
+            road_sec_dens = road_sec_len / window_len
+            road_4wd_dens = road_4wd_len / window_len
 
         curs.execute(f'UPDATE IGOAttributes SET Road_len = {road_len} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET Rail_len = {rail_len} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET Canal_len = {canal_len} WHERE IGOID = {igoid}')
-        curs.execute(f'UPDATE IGOAttributes SET Roadx_ct = {roadx_ct} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET RoadX_ct = {roadx_ct} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET DivPts_ct = {divpts_ct} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET Road_prim_len = {road_prim_len} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET Road_sec_len = {road_sec_len} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET Road_4wd_len = {road_4wd_len} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET Road_dens = {road_dens} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET Rail_dens = {rail_dens} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET Canal_dens = {canal_dens} WHERE IGOID = {igoid}')
-        curs.execute(f'UPDATE IGOAttributes SET Roadx_dens = {roadx_dens} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET RoadX_dens = {roadx_dens} WHERE IGOID = {igoid}')
         curs.execute(f'UPDATE IGOAttributes SET DivPts_dens = {divpts_dens} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET Road_prim_dens = {road_prim_dens} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET Road_sec_dens = {road_sec_dens} WHERE IGOID = {igoid}')
+        curs.execute(f'UPDATE IGOAttributes SET Road_4wd_dens = {road_4wd_dens} WHERE IGOID = {igoid}')
 
     conn.commit()
     conn.close()

packages/anthro/database/anthro_schema.sql

@@ -58,7 +58,8 @@ CREATE TABLE Watersheds (
 
 CREATE TABLE IGOAttributes (
     IGOID INTEGER PRIMARY KEY NOT NULL,
-    LevelPathI REAL,
+    FCode INTEGER,
+    level_path REAL,
     seg_distance REAL,
     stream_size INTEGER,
     LUI REAL,
@@ -71,11 +72,18 @@ CREATE TABLE IGOAttributes (
     RoadX_ct INTEGER,
     RoadX_dens REAL,
     DivPts_ct INTEGER,
-    DivPts_dens REAL);
+    DivPts_dens REAL,
+    Road_prim_len REAL,
+    Road_prim_dens REAL,
+    Road_sec_len REAL,
+    Road_sec_dens REAL,
+    Road_4wd_len REAL,
+    Road_4wd_dens REAL);
 
 CREATE TABLE DGOAttributes (
     DGOID INTEGER PRIMARY KEY NOT NULL,
-    LevelPathI REAL,
+    FCode INTEGER,
+    level_path REAL,
     seg_distance REAL,
     centerline_length REAL,
     segment_area REAL,
@@ -84,7 +92,10 @@ CREATE TABLE DGOAttributes (
     Rail_len REAL,
     Canal_len REAL,
     RoadX_ct INTEGER,
-    DivPts_ct INTEGER);
+    DivPts_ct INTEGER,
+    Road_prim_len REAL,
+    Road_sec_len REAL,
+    Road_4wd_len REAL);
 
 CREATE TABLE ReachAttributes (
     ReachID INTEGER PRIMARY KEY NOT NULL,

packages/anthro/database/data/LandUses.csv

@@ -53,15 +53,15 @@ LandUseID,Name,Intensity
 52,Depressional Wetland,0
 53,Desert Scrub,0
 54,Developed-High Intensity,1
-55,Developed-Low Intensity,1
+55,Developed-Low Intensity,0.66
 56,Developed-Medium Intensity,1
 57,Developed-Open Space,0.33
 58,Developed-Roads,1
-59,Developed-Upland Deciduous Forest,1
-60,Developed-Upland Evergreen Forest,1
-61,Developed-Upland Herbaceous,1
-62,Developed-Upland Mixed Forest,1
-63,Developed-Upland Shrubland,1
+59,Developed-Upland Deciduous Forest,0.66
+60,Developed-Upland Evergreen Forest,0.66
+61,Developed-Upland Herbaceous,0.66
+62,Developed-Upland Mixed Forest,0.66
+63,Developed-Upland Shrubland,0.66
 64,Douglas-fir Forest and Woodland,0
 65,Douglas-fir-Grand Fir-White Fir Forest and Woodland,0
 66,Douglas-fir-Ponderosa Pine-Lodgepole Pine Forest and Woodland,0

scripts/automation/FargateAnthro.sh

@@ -82,7 +82,7 @@ try() {
     $RS_CONTEXT_DIR/topography/dem_hillshade.tif \
     $VBET_DIR/outputs/vbet.gpkg/vbet_igos \
     $VBET_DIR/intermediates/vbet_intermediates.gpkg/vbet_dgos \
-    $RS_CONTEXT_DIR/hydrology/hydro_derivatives.gpkg/network_intersected_300m \
+    $RS_CONTEXT_DIR/hydrology/hydro_derivatives.gpkg/network_intersected \
     $VBET_DIR/outputs/vbet.gpkg/vbet_full \
     $RS_CONTEXT_DIR/ownership/ownership.shp \
     $RS_CONTEXT_DIR/transportation/canals.shp \