code cleanup

packages/anthro/anthro/anthro.py

@@ -23,8 +23,7 @@
 from rscommons.augment_lyr_meta import augment_layermeta, add_layer_descriptions, raster_resolution_meta
 
 from anthro.utils.conflict_attributes import conflict_attributes
-# from anthro.igo_infrastructure import infrastructure_attributes
-from anthro.utils.igo_infrastructure2 import infrastructure_attributes
+from anthro.utils.igo_infrastructure import infrastructure_attributes
 from anthro.utils.igo_vegetation import igo_vegetation
 from anthro.utils.igo_land_use import calculate_land_use
 from anthro.utils.lui_raster import lui_raster
@@ -168,8 +167,6 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
     # Execute the SQL to create the lookup tables in the output geopackage
     create_database(huc, outputs_gpkg_path, db_metadata, cfg.OUTPUT_EPSG, os.path.join(os.path.abspath(os.path.dirname(__file__)), '..', 'database', 'anthro_schema.sql'))
 
-    # project.add_metadata_simple(db_metadata)
-
     igo_geom_path = os.path.join(outputs_gpkg_path, LayerTypes['OUTPUTS'].sub_layers['ANTHRO_GEOM_POINTS'].rel_path)
     line_geom_path = os.path.join(outputs_gpkg_path, LayerTypes['OUTPUTS'].sub_layers['ANTHRO_GEOM_LINES'].rel_path)
     dgo_geom_path = os.path.join(outputs_gpkg_path, LayerTypes['OUTPUTS'].sub_layers['ANTHRO_GEOM_DGOS'].rel_path)
@@ -233,13 +230,7 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
          RSMeta('Very Large Search Window', str(distancein['3']), RSMetaTypes.INT, locked=True),
          RSMeta('Huge Search Window', str(distancein['4']), RSMetaTypes.INT, locked=True)])
 
-    # get moving window for each igo
-    # windows = get_moving_windows(igo_geom_path, input_layers['DGO'], levelpathsin, distancein)
-    # with SQLiteCon(outputs_gpkg_path) as database:
-    #     for fid, windowvals in windows.items():
-    #         database.curs.execute(f'UPDATE IGOAttributes SET window_area = {windowvals[2]} WHERE IGOID = {fid}')
-    #         database.curs.execute(f'UPDATE IGOAttributes SET window_length = {windowvals[1]} WHERE IGOID = {fid}')
-    #     database.conn.commit()
+    # associate DGO IDs with IGO IDs for moving windows
     windows = moving_window_dgo_ids(igo_geom_path, input_layers['DGO'], levelpathsin, distancein)
 
     # calculate conflict attributes for reaches
@@ -252,8 +243,6 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
     infrastructure_attributes(windows, input_layers['ROADS'], input_layers['RAILS'], input_layers['CANALS'],
                               crossings, diversions, outputs_gpkg_path)
 
-    # add layers to project
-
     # get land use attributes for reaches
     vegetation_summary(outputs_gpkg_path, input_layers['DGO'], existing_veg)
     land_use(outputs_gpkg_path)
@@ -262,22 +251,11 @@ def anthro_context(huc: int, existing_veg: Path, hillshade: Path, igo: Path, dgo
     calculate_land_use(outputs_gpkg_path, windows)
     lui_raster(existing_veg, outputs_gpkg_path, os.path.join(os.path.dirname(intermediates_gpkg_path), 'lui.tif'))
     # add lui raster to project
-    # project.add_dataset(proj_nodes['Intermediates'], os.path.join(os.path.dirname(intermediates_gpkg_path), 'lui.tif'), LayerTypes['LUI'], 'Raster')
     lui_node, lui_ras = project.add_project_raster(proj_nodes['Intermediates'], LayerTypes['LUI'])
     raster_resolution_meta(project, lui_ras, lui_node)
 
     ellapsed_time = time.time() - start_time
 
-    # intermediate_layers = [
-    #     ['diversions', 'Stream Diversion Points'],
-    #     ['private_land', 'Private Land'],
-    #     ['rail_valleybottom', 'Railroads in Valley Bottom'],
-    #     ['road_crossings', 'Road Stream Crossings'],
-    #     ['road_valleybottom', 'Roads in Valley Bottom']
-    # ]
-    # for lyrname in intermediate_layers:
-    #     LayerTypes['INTERMEDIATES'].add_sub_layer(f'{str.upper(lyrname[0])}', RSLayer(f'{lyrname[1]}', f'{str.upper(lyrname[0])}', 'Vector', f'{lyrname[0]}'))
-
     project.add_project_geopackage(proj_nodes['Intermediates'], LayerTypes['INTERMEDIATES'])
     project.add_metadata([
         RSMeta("ProcTimeS", "{:.2f}".format(ellapsed_time), RSMetaTypes.HIDDEN, locked=True),

packages/anthro/anthro/utils/conflict_attributes.py

@@ -63,8 +63,8 @@ def conflict_attributes(
 
 def calc_conflict_attributes(flowlines_path, valley_bottom, roads, rail, canals, ownership, buffer_distance_metres, cell_size_meters, epsg, canal_codes, intermediates_gpkg_path):
 
-    log = Logger('Conflict')
-    log.info('Calculating conflict attributes')
+    log = Logger('Infrastructure')
+    log.info('Calculating distance from channel to infrastructure')
 
     # Create union of all reaches and another of the reaches without any canals
     reach_union = get_geometry_unary_union(flowlines_path)
@@ -90,31 +90,13 @@ def calc_conflict_attributes(flowlines_path, valley_bottom, roads, rail, canals,
 
     # Buffer all reaches (being careful to use the units of the Shapefile)
     reaches = load_geometries(flowlines_path, epsg=epsg)
-    # dgo_geoms = load_geometries(dgos_path, epsg=epsg)
+
     geopackage_path = os.path.dirname(flowlines_path)
     with get_shp_or_gpkg(ownership) as lyr:
         buffer_distance = lyr.rough_convert_metres_to_vector_units(buffer_distance_metres)
         cell_size = lyr.rough_convert_metres_to_vector_units(cell_size_meters)
 
-    # st = datetime.datetime.now()
-    polygons = {}
-    # for reach_id, polyline in reaches.items():
-    #     log.info(f'finding IGOs that intersect network segment {reach_id}')
-    #     polys = []
-    #     for dgo_id, polygon in dgo_geoms.items():
-    #         if polygon.intersects(polyline):
-    #             polys.append(polygon)
-    #     if len(polys) > 1:
-    #         poly = unary_union(polys)
-    #     elif len(polys) == 1:
-    #         poly = polys[0]
-    #     else:  # if there are no polygons that intersect network?
-    #         poly = polyline.buffer(buffer_distance)
-    #     polygons[reach_id] = poly
-
     polygons = {reach_id: polyline.buffer(buffer_distance) for reach_id, polyline in reaches.items()}
-    # end = datetime.datetime.now()
-    # print(f'finding intersecting dgos took {end-st}')
 
     results = {}
     tmp_folder = os.path.join(os.path.dirname(intermediates_gpkg_path), 'tmp_conflict')
@@ -140,7 +122,7 @@ def calc_conflict_attributes(flowlines_path, valley_bottom, roads, rail, canals,
     if len(reaches) > 0:
         admin_agency(geopackage_path, reaches, ownership, results)
 
-    log.info('Conflict attribute calculation complete')
+    log.info('Infrastructure attribute calculation complete')
 
     # Cleanup temporary feature classes
     safe_remove_dir(tmp_folder)

packages/anthro/anthro/utils/igo_infrastructure.py

@@ -1,18 +1,23 @@
 """Functions to attribute IGO points with attributes related to the presence of infrastructure within the riverscape.
-
 Jordan Gilbert
-
 Dec 2022
 """
+import os
 import sqlite3
-from rscommons import get_shp_or_gpkg
+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
 from rscommons.vector_ops import get_geometry_unary_union
 
 
-def infrastructure_attributes(igo: str, windows: str, road: str, rail: str, canal: str, crossings: str, diversions: str,
+def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, crossings: str, diversions: str,
                               out_gpkg_path: str):
 
+    log = Logger('IGO Infrastructure Attributes')
+    log.info('Adding attributes for infrastructure density to IGOs')
+
     in_data = {
         road: 'Road',
         rail: 'Rail',
@@ -22,66 +27,123 @@ def infrastructure_attributes(igo: str, windows: str, road: str, rail: str, cana
     }
 
     # get epsg_proj
-    with get_shp_or_gpkg(igo) as inref:
+    with get_shp_or_gpkg(road) as inref:
         ftr = inref.ogr_layer.GetFeature(1)
-        long = ftr.GetGeometryRef().GetEnvelope()[0]
-        epsg_proj = get_utm_zone_epsg(long)
+        if ftr is not None:
+            long = ftr.GetGeometryRef().GetEnvelope()[0]
+            epsg_proj = get_utm_zone_epsg(long)
+        else:
+            epsg_proj = None
+
+    if not epsg_proj:
+        with get_shp_or_gpkg(road) 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)
+            else:
+                epsg_proj = None
+
+    if not epsg_proj:
+        with get_shp_or_gpkg(canal) 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)
+            else:
+                epsg_proj = None
+
+    if not epsg_proj:
+        log.info('No infrastructure datasets, skipping attributes')
+
+        return
 
     conn = sqlite3.connect(out_gpkg_path)
+    curs = conn.cursor()
 
-    for dataset, label in in_data.items():
-        # ds = get_geometry_unary_union(dataset)
-        # sref, transform = ds.get_transform_from_epsg(ds.spatial_ref, epsg_proj)
-        print(f'Summarizing metrics for dataset: {dataset}')
+    with get_shp_or_gpkg(road) as reflyr:
+        sref, transform = reflyr.get_transform_from_epsg(reflyr.spatial_ref, epsg_proj)
 
-        counter = 1
-        for igoid, window in windows.items():
-            # print(f'summarizing on igo {counter} of {len(windows)} for dataset {label}')
+    gpkg_driver = ogr.GetDriverByName('GPKG')
+    with GeopackageLayer(out_gpkg_path, 'DGOGeometry') as dgo_lyr:
+        for dataset, label in in_data.items():
+            log.info(f'Calculating metrics for dataset: {label}')
+            dsrc = gpkg_driver.Open(os.path.dirname(dataset))
+            if dsrc.GetLayer(os.path.basename(dataset)) is None:
+                continue
+            ds = get_geometry_unary_union(dataset)
+            if ds is None:
+                log.info(f'Skipping dataset {label} because it contains no features')
+                continue
 
-            ftr_length = 0.0
-            ftr_count = 0
+            for dgo_ftr, *_ in dgo_lyr.iterate_features(f'summarizing dataset {label} metrics on DGOs'):
+                dgoid = dgo_ftr.GetFID()
 
-            if igoid in [1110, 1442]:
-                print(igoid)
+                dgo_ogr = dgo_ftr.GetGeometryRef()
+                dgo_geom = VectorBase.ogr2shapely(dgo_ogr)
 
-            with get_shp_or_gpkg(dataset) as layer:
-                sref, transform = layer.get_transform_from_epsg(layer.spatial_ref, epsg_proj)
-                for geom in window[0].geoms:
-                    layer.ogr_layer.SetSpatialFilter(VectorBase.shapely2ogr(geom))
-                    if layer.ogr_layer.GetFeatureCount() == 0:
-                        continue
-                    else:
-                        for feature in layer.ogr_layer:
-                            lyr_cl = geom.intersection(VectorBase.ogr2shapely(feature))
-
-                            # leave null if layer is empty?
-                            if lyr_cl.is_empty is True:
-                                continue
-                            else:
-                                if lyr_cl.type in ['MultiLineString', 'LineString']:
-                                    ogrlyr = VectorBase.shapely2ogr(lyr_cl)
-                                    lyr_clipped = VectorBase.ogr2shapely(ogrlyr, transform=transform)
-                                    ftr_length += lyr_clipped.length
-                                if lyr_cl.type in ['MultiPoint']:
-                                    ftr_count += len(lyr_cl.geoms)
-                                if lyr_cl.type in ['Point']:
-                                    if lyr_cl.is_empty is False:
-                                        ftr_count += 1
-
-                if layer.ogr_geom_type in layer.LINE_TYPES:
+                lyr_cl = dgo_geom.intersection(ds)
+
+                if lyr_cl.type in ['MultiLineString', 'LineString']:
                     lb1 = label + '_len'
-                    lb2 = label + '_dens'
-                    conn.execute(f'UPDATE IGOAttributes SET {lb1} = {ftr_length} WHERE IGOID = {igoid}')
-                    if window[2] != 0.0:
-                        conn.execute(f'UPDATE IGOAttributes SET {lb2} = {ftr_length / window[2]} WHERE IGOID = {igoid}')
+                    if lyr_cl.is_empty is True:
+                        curs.execute(f'UPDATE DGOAttributes SET {lb1} = 0 WHERE DGOID = {dgoid}')
+                    else:
+                        ogrlyr = VectorBase.shapely2ogr(lyr_cl)
+                        lyr_clipped = VectorBase.ogr2shapely(ogrlyr, transform=transform)
+                        curs.execute(f'UPDATE DGOAttributes SET {lb1} = {lyr_clipped.length} WHERE DGOID = {dgoid}')
 
-                if layer.ogr_geom_type in layer.POINT_TYPES:
+                if lyr_cl.type in ['MultiPoint']:
                     lb1 = label + '_ct'
-                    lb2 = label + '_dens'
-                    conn.execute(f'UPDATE IGOAttributes SET {lb1} = {ftr_count} WHERE IGOID = {igoid}')
-                    if window[2] != 0.0:
-                        conn.execute(f'UPDATE IGOAttributes SET {lb2} = {ftr_count / window[2]} WHERE IGOID = {igoid}')
-
-                conn.commit()
+                    if lyr_cl.is_empty is True:
+                        curs.execute(f'UPDATE DGOAttributes SET {lb1} = 0 WHERE DGOID = {dgoid}')
+                    else:
+                        curs.execute(f'UPDATE DGOAttributes SET {lb1} = {len(lyr_cl.geoms)} WHERE DGOID = {dgoid}')
 
-                counter += 1
+                if lyr_cl.type in ['Point']:
+                    lb1 = label + '_ct'
+                    if lyr_cl.is_empty is True:
+                        curs.execute(f'UPDATE DGOAttributes SET {lb1} = 0 WHERE DGOID = {dgoid}')
+                    else:
+                        curs.execute(f'UPDATE DGOAttributes SET {lb1} = 1 WHERE DGOID = {dgoid}')
+
+    fields = ['Road_len', 'Rail_len', 'Canal_len', 'RoadX_ct', 'DivPts_ct']
+    for field in fields:
+        curs.execute(f'UPDATE DGOAttributes SET {field} = 0 WHERE {field} IS NULL')
+
+    # summarize metrics from DGOs to IGOs using moving windows
+    for igoid, dgoids in windows.items():
+        road_len = 0
+        rail_len = 0
+        canal_len = 0
+        roadx_ct = 0
+        divpts_ct = 0
+        window_area = 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_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
+
+        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 DivPts_ct = {divpts_ct} 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 DivPts_dens = {divpts_dens} WHERE IGOID = {igoid}')
+
+    conn.commit()
+    conn.close()