Lee Hachadoorian wrote:
> I'm looking through the PostGIS reference, and I can't seem to find a
> way to take a geometry of polygons and turn it into lines. What I'm
> looking for is something like the ArcGIS Feature to Line geoprocessor,
> which will create a line shapefile where each feature is an arc
> representing the boundary between neighboring polygons with a field
> indicating the ids of the polygon on either side.
>
> Functions like ST_MakeLine require point geometries, and I don't see
> anything else that seems to be what I'm looking for. Any ideas would
> be welcome.
>
> Thanks,
> Lee Hachadoorian
> PhD Student in Geography
> Program in Earth & Environmental Sciences
> CUNY Graduate Center
A custom function is the fastest and most efficient way to do this, as others
have indicated,
but it's often fun and instructive to try to implement a task in a single
statement.
Here's my attempt at the "ArcGIS Feature to Line geoprocessor" mentioned.
It needs to be faster and smarter, but seems to work. If anyone has
suggestions for this or even a different
approach, please chime in. I'm not always good at seeing other paths once I
get a ways into the woods.
Thanks.
-Eric
/*
turn_poly_into_lines.sql
I used parcel polygons in my experiment.
Replace "polytable" with your poly table. Replace gid value with your
poly's gid.
*/
select p4.poly_id,p4.adjacent_poly_id,p4.pos_order,p4.geom
from
(
/* Get edges shared with other polys */
select t1.gid as poly_id,t2.gid as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as pos_order,
ST_linemerge(ST_intersection(ST_linemerge(ST_boundary(t1.geom)),ST_linemerge(ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
union
/*
Create edges not shared with other polys. Uses -1 for pseudo external
universe polygon.
Produces extra line necessitating the self-join and the ST_covers filter at
the end.
Need help, more work on this.
*/
select p1.poly_id, p1.adjacent_poly_id, p1.position as
pos_order,ST_line_substring(p3.geom,p1.position,p2.position) as geom
from
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as position,
ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and
ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
not in
(
select
ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p1,
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as position,
ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and
ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
not in
(
select
ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p2,
(
select ST_linemerge(ST_boundary(geom)) as geom
from polytable
where gid = 309
)
as p3
where p1.position < p2.position
)
as p4,
/* Repeat for needed self-join */
(
/* Get edges shared with other polys */
select t1.gid as poly_id,t2.gid as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as pos_order,
ST_linemerge(ST_intersection(ST_linemerge(ST_boundary(t1.geom)),ST_linemerge(ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
union
/*
Create edges not shared with other polys. Uses -1 for pseudo external
universe polygon.
Produces extra line necessitating the self-join and the ST_covers filter at
the end.
Need help, more work on this.
*/
select p1.poly_id, p1.adjacent_poly_id, p1.position as pos_order,
ST_line_substring(p3.geom,p1.position,p2.position) as geom
from
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as position,
ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and
ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
not in
(
select
ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p1,
(
select t1.gid as poly_id,-1::integer as adjacent_poly_id,
ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
as position,
ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
and
ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
not in
(
select
ST_endpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))
as geom
from polytable t1, polytable t2
where t1.gid = 309
and ST_touches(t1.geom, t2.geom)
)
)
as p2,
(
select ST_linemerge(ST_boundary(geom)) as geom
from polytable
where gid = 309
)
as p3
where p1.position < p2.position
)
as p5
where ST_covers(p4.geom,p5.geom)
group by p4.poly_id,p4.adjacent_poly_id,p4.pos_order,p4.geom
having count(ST_covers(p4.geom,p5.geom)) = 1
order by p4.pos_order
---------
Eric Randall
County of Erie
140 W 6th St
Room 119
Erie, PA 16501
ph. 814-451-6063
fx. 814-451-7000
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
