Skip to content

Commit

Permalink
density by road type
Browse files Browse the repository at this point in the history
  • Loading branch information
@jtgilbert @philipbaileynar
jtgilbert authored and philipbaileynar committed Feb 2, 2024
1 parent 0db0552 commit e89368a
Show file tree
Hide file tree
Showing 2 changed files with 70 additions and 7 deletions.
64 changes: 59 additions & 5 deletions packages/anthro/anthro/utils/igo_infrastructure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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

Expand All @@ -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))
Expand Down Expand Up @@ -108,6 +110,37 @@ 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()

Expand All @@ -117,6 +150,9 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
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 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():
Expand All @@ -125,6 +161,9 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
canal_len = 0
roadx_ct = 0
divpts_ct = 0
road_prim_len = 0
road_sec_len = 0
road_4wd_len = 0
window_len = 0

for dgoid in dgoids:
Expand All @@ -133,6 +172,9 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
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]
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_len == 0:
Expand All @@ -141,23 +183,35 @@ def infrastructure_attributes(windows: str, road: str, rail: str, canal: str, cr
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_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 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 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()
13 changes: 11 additions & 2 deletions packages/anthro/database/anthro_schema.sql
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,13 @@ 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,
Expand All @@ -86,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,
Expand Down

0 comments on commit e89368a

Please sign in to comment.