Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Sketch for topological simplification of national hydrology dataset
'''
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'))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.