Sketch for topological simplification of national hydrology dataset

snip-lines.py

'''

ogr2ogr -overwrite -t_srs '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs' PG:"user=nhd dbname=nhd" NHDFlowline.shp

create table nhdflowline_snapped ( ogc_fid integer primary key );

select addgeometrycolumn('', 'nhdflowline_snapped', 'geometry', 900914, 'LINESTRING', 2);
create index nhdflowline_snapped_geom_idx on nhdflowline_snapped using gist (geometry);

insert into nhdflowline_snapped
(ogc_fid, geometry)
select ogc_fid, snaptogrid(force_2d(wkb_geometry), 1)
from nhdflowline



ogr2ogr -overwrite -t_srs '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs' PG:"user=nhd dbname=nhd" NHDArea.shp

create table nhdarea_snapped ( ogc_fid integer primary key );

select addgeometrycolumn('', 'nhdarea_snapped', 'geometry', 900914, 'POLYGON', 2);
create index nhdarea_snapped_geom_idx on nhdarea_snapped using gist (geometry);

insert into nhdarea_snapped
(ogc_fid, geometry)
select ogc_fid, snaptogrid(force_2d(wkb_geometry), 1)
from nhdarea

'''

import json
from itertools import groupby
from psycopg2 import connect
from shapely import wkb
from shapely.geometry import Polygon

if __name__ == '__main__':
    
    db = connect(user='nhd', database='nhd').cursor()
    
    #
    # Look at each NHD Area in turn by ID.
    #
    q = '''SELECT ogc_fid
           FROM nhdarea_snapped
           ORDER BY ogc_fid'''
    
    db.execute(q)
    
    features = []
    
    for (afid, ) in db.fetchall():
    
        #
        # Chop a single area into portions of boundary linestring and
        # locations (0 - 1 inclusive) along those linestrings showing
        # breakpoints where smaller streams touch the surface.
        #
        q = '''SELECT (dump(boundary(a.geometry))).path AS path,
                      asbinary((dump(boundary(a.geometry))).geom) AS geom,
                      line_locate_point((dump(boundary(a.geometry))).geom, intersection(a.geometry, fl.geometry)) AS loc
               FROM nhdarea_snapped AS a,
                    nhdflowline_snapped AS fl
               WHERE a.ogc_fid = %s
                 AND a.geometry && fl.geometry
                 AND touches(a.geometry, fl.geometry)
               ORDER BY path, loc'''
        
        db.execute(q, (afid, ))
        
        rings = []
        
        #
        # Iterate, grouping by part path to form a complete ring.
        #
        for (path, parts) in groupby(db.fetchall(), lambda row: row[0]):
        
            parts = list(parts)
            shape = parts[0][1]
            locs = [part[2] for part in parts]
            
            if locs[0] != 0.0:
                locs.insert(0, 0.0)
        
            if locs[-1] != 1.0:
                locs.append(1.0)
        
            ring = []
            
            #
            # Simplify each portion of the overall linestring in turn.
            #
            for (loc1, loc2) in zip(locs, locs[1:]):
            
                q = '''SELECT asbinary(simplify(line_substring(geomfromwkb(%s), %s, %s), 2.5))'''
                
                db.execute(q, (shape, loc1, loc2))
                row = db.fetchone()
                
                ring.extend(wkb.loads(str(row[0])).coords)
                
            rings.append(ring)
        
        if len(rings) == 0:
            continue
        
        feature = dict(type='Feature', id=afid, properties=dict(ogc_fid=afid))
        feature['geometry'] = Polygon(rings[0], rings[1:]).__geo_interface__
        
        features.append(feature)

    json.dump(dict(type='FeatureCollection', features=features), open('out.json', 'w'))