[postgis-users] Turn a Polygon into Lines

Randall, Eric ERandall at eriecountygov.org
Sun Apr 6 08:47:05 PDT 2008


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



This is a correction to my March 31 post.   I should have have tested a little more and not been in such a hurry I guess.
I neglected to account for the cases where a non-shared edge needs to end at the beginning.  Adding a case statement
seems to have handled that.   I'd like to find a way to reference the derived table (p4) multiple times to avoid
repeating the query  (p5).  The Common Table Expression WITH statement should do this I think but it doesn't seem to be 
supported in my version of PostgreSQL (8.2).  Can anyone give me some ideas?  Thanks.

-Eric



/* 

    example_turn_poly_into_lines.sql
   
   replace "polytable" with your poly table.
   replace gid value (309) 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(s) where there is more than one non-shared edge,
   necessitating the self-join and ST_covers filter at the end.
   Need help, more work on this, ...suggestions?
*/

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, 
case ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))) 
when 0 then 1
else ST_line_locate_point(ST_linemerge(boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
end
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(s) where there is more than one non-shared edge,
   necessitating the self-join and ST_covers filter at the end.
   Need help, more work on this, ...suggestions?
*/

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, 
case ST_line_locate_point(ST_linemerge(ST_boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom))))) 
when 0 then 1
else ST_line_locate_point(ST_linemerge(boundary(t1.geom)),ST_startpoint(ST_linemerge(ST_intersection(ST_boundary(t1.geom),ST_boundary(t2.geom)))))
end
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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20080406/5b5d2a77/attachment.html>


More information about the postgis-users mailing list