[postgis-users] Getting dimensions of polygon sides

Simon Greener simon at spatialdbadvisor.com
Thu Jun 12 17:43:38 PDT 2008


I have created some functions in Oracle that "vectorises" a linear or polygon geometry that I use in centroid determination and also for cartographic purposes (eg the display of a bearing and distance for each side of a land parcel).

I am no PostGIS expert and enjoy lurking on this list as I learn a lot of good things from the experts that contribute to it. But the vectorising of a geometry was a nice exercise to see if I could port my Oracle code to PostGIS and learn something along the way. The following code was created months ago and is posted now in case anyone finds it useful. (I would also appreciate pointers from the experts here on how to handled geometries with curved elements.)

CREATE TYPE coordtype AS
    (x double precision,
     y double precision,
     z double precision,
     m double precision);

CREATE TYPE vectortype AS
    (startcoord coordtype,
     endcoord coordtype);

-------------------------------------------------------------
DROP              FUNCTION GetVectorViaSQL(p_geometry geometry);
CREATE OR REPLACE FUNCTION GetVectorViaSQL(p_geometry geometry)
   RETURNS SETOF VectorType STABLE AS $$
DECLARE
     v_GeometryType   varchar(1000);
     v_rec            RECORD;
     v_vector         VectorType;
     v_start          CoordType;
     v_end            CoordType;
     c_linestring CURSOR ( p_geom geometry )
     IS
       SELECT X(sp) as sx,Y(sp) as sy,Z(sp) as sz,M(sp) as sm,X(ep) as ex,Y(ep) as ey,Z(ep) as ez,M(ep) as em
        FROM ( SELECT pointn(p_geom, generate_series(1, npoints(p_geom)-1)) as sp,
                      pointn(p_geom, generate_series(2, npoints(p_geom)  )) as ep
             ) AS foo;
     c_multilinestring CURSOR ( p_geom geometry )
     IS
       SELECT X(sp),Y(sp),Z(sp),M(sp),X(ep),Y(ep),Z(ep),M(ep)
         FROM ( SELECT pointn(b.geom, generate_series(1, npoints(b.geom)-1)) as sp,
                       pointn(b.geom, generate_series(2, npoints(b.geom)  )) as ep
                  FROM (select geometryn(p_geom,generate_series(1,numgeometries(p_geom))) as geom) as b
              ) AS foo;
     c_polygon CURSOR ( p_geom geometry )
     IS
       SELECT X(sp),Y(sp),Z(sp),M(sp),X(ep),Y(ep),Z(ep),M(ep)
         FROM ( SELECT pointn(a.geom, generate_series(1, npoints(a.geom)-1)) as sp,
                       pointn(a.geom, generate_series(2, npoints(a.geom)  )) as ep
                  FROM ( SELECT exteriorring(p_geom) as geom
                         UNION ALL
                         SELECT interiorringn(p_geom,generate_series(1,numinteriorrings(p_geom))) as geom
                       ) a
              ) b;
     c_multipolygon CURSOR ( p_geom geometry )
     IS
       SELECT X(sp),Y(sp),Z(sp),M(sp),X(ep),Y(ep),Z(ep),M(ep)
         FROM ( SELECT pointn(a.geom, generate_series(1, npoints(a.geom)-1)) as sp,
                       pointn(a.geom, generate_series(2, npoints(a.geom)  )) as ep
                  FROM ( SELECT exteriorring(b.geom) as geom
                           FROM (select geometryn(p_geom,generate_series(1,numgeometries(p_geom))) as geom ) b
                         UNION ALL
                         SELECT interiorringn(c.geom,generate_series(1,numinteriorrings(c.geom))) as geom
                           FROM (select geometryn(p_geom,generate_series(1,numgeometries(p_geom))) as geom ) c
                       ) a
                     ) b;
     c_refcursor refcursor;
Begin
     If ( p_geometry is NULL ) Then
       return;
     End If;
     v_GeometryType := ST_GeometryType(p_geometry);
     If    ( v_GeometryType in ('ST_Point','ST_MultiPoint') ) Then
       return;
     End If;
     IF    ( v_GeometryType = 'ST_MultiLineString' ) THEN
       OPEN c_multilinestring(p_geometry);
       c_refcursor := c_multilinestring;
     ELSIF ( v_GeometryType = 'ST_LineString'      ) THEN
       OPEN c_linestring(p_geometry);
       c_refcursor := c_linestring;
     ELSIF ( v_GeometryType = 'ST_MultiPolygon'    ) THEN
       OPEN c_multipolygon(p_geometry);
       c_refcursor := c_multipolygon;
     ELSIF ( v_GeometryType = 'ST_Polygon'         ) THEN
       OPEN c_polygon(p_geometry);
       c_refcursor := c_polygon;
     ELSIF ( v_GeometryType = 'ST_Geometry'        ) THEN
       -- COuld be anything... use native PostGIS function
       v_GeometryType = GeometryType(p_geometry);
       IF ( v_geometryType = 'GEOMETRYCOLLECTION' ) THEN
          FOR v_geom IN 1..ST_NumGeometries(p_geometry) LOOP
             FOR v_rec IN SELECT * FROM GetVectorViaSQL(ST_GeometryN(p_geometry,v_geom)) LOOP
    	        RETURN NEXT v_rec;
    	    END LOOP;
          END LOOP;
       ELSE
          -- Probably CURVED something...
          RETURN;
       END IF;
     END IF;
     IF ( v_geometryType NOT IN ('ST_Geometry','GEOMETRYCOLLECTION') ) THEN
       LOOP
         FETCH c_refcursor INTO
               v_start.x, v_start.y, v_start.z, v_start.m,
               v_end.x,   v_end.y,   v_end.z,   v_end.m;
         v_vector.startcoord := v_start;
         v_vector.endcoord   := v_end;
         EXIT WHEN NOT FOUND;
         RETURN NEXT v_vector;
       END LOOP;
       CLOSE c_refcursor;
     END IF;
end;
$$ LANGUAGE plpgsql;

select * from GetVectorViaSQL('GEOMETRYCOLLECTION(POINT(2 3 4),LINESTRING(2 3 4,3 4 5))'::geometry);
select * from GetVectorViaSQL('LINESTRING(0 0, 1 1, 2 2, 3 3)'::geometry) as geom;
select * from GetVectorViaSQL('MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))'::geometry) As GV;
select * from GetVectorViaSQL('POLYGON((326454.7 5455793.7,326621.3 5455813.7,326455.4 5455796.6,326454.7 5455793.7))'::geometry);
select * from GetVectorViaSQL('MULTIPOLYGON(((326454.7 5455793.7,326621.3 5455813.7,326455.4 5455796.6,326454.7 5455793.7)),((326771.6 5455831.6,326924.1 5455849.9,326901.9 5455874.2,326900.7 5455875.8,326888.9 5455867.3,326866 5455853.1,326862 5455851.2,326847.4 5455845.8,326827.7 5455841.2,326771.6 5455831.6)))'::geometry);

regards
Simon
On Fri, 13 Jun 2008 04:03:15 +1000, Obe, Regina <robe.dnd at cityofboston.gov> wrote:

> Bruce,
> Yes that is what I was thinking.  I was going to create another geometry field which is a point which represents the exact location I want the annotation placed.  I was using quantum for the below snapshots and I think it just places it at the centroid of the line.  I looked at it in ArcPad which is what will be used in the field and it seems to do more or less the same too although I think it auto angles so I may go with instead of centroid - just another line translated 90 degrees some x amount to the right to maintain the angularness of the text.
>Now if only I can remember my highschool trig.   What you are saying makes perfect sense though so I'm sure I will be able to figure out the functional steps with a little bit of thinking.
>Thanks,
> Regina
>
> ________________________________
>
> From: postgis-users-bounces at postgis.refractions.net on behalf of Bruce Rindahl
> Sent: Thu 6/12/2008 1:03 PM
> To: PostGIS Users Discussion
> Subject: Re: [postgis-users] Getting dimensions of polygon sides
>
>
>
> Obe, Regina wrote:
>> Anyrate my new problem:  This whole exercise just made me realize they
>> must have put these polygons together with the constituent linework of
>> the edges.  Because if I look at the abutting parcels those extra dims
>> represent sides of those smaller parcels.  My simplify solution works
>> great if you are just looking at one parcel, but when you look at
>> abutters - it shall we say, looks crowded.  So my thought is to put the
>> position of my length annotation inside each parcel it represents.
>>
>> I'm lost how to do that.  So basically jimmy the ST_Centroid(the_geom)
>> of my above such that it always sits inside the parcel boundary.
>>
>> Thanks,
>> Regina
>>
> I am not sure how you are placing the annotation, but this could be
> interesting....
>
> Looking at the aftersimplify.png, you have the length of each line
> segment (and the segment itself).  The center point of each segment is
> easy.  Now you need to define a line segment offset by the text height
> of your annotation and translate it INSIDE the polygon.  The you draw
> (label?, annotate? ....) the text along that line.  If you first update
> the polygon with forceRHR() then INSIDE is always to the right of the
> line segment.  You can easily get the azmith of the linesegment from the
> start and end points, add 90 degrees and compute the x and y values for
> the translate() function.  This should be basic high school trig - ;-) .
>
> Is this what you are thinking of?
>
> Bruce
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
>
>
> -----------------------------------------
> The substance of this message, including any attachments, may be
> confidential, legally privileged and/or exempt from disclosure
> pursuant to Massachusetts law. It is intended
> solely for the addressee. If you received this in error, please
> contact the sender and delete the material from any computer.
>



-- 
SpatialDB Advice and Design, Solutions Architecture and Programming,
Oracle Database 10g Administrator Certified Associate; Oracle Database 10g SQL Certified Professional
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE, Manifold GIS, Radius Topology and Studio Specialist.
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia.
Website: www.spatialdbadvisor.com
   Email: simon at spatialdbadvisor.com
   Voice: +613 9016 3910
Mobile: +61 418 396391
Skype: sggreener
Longitude: 147.20515 (147° 12' 18" E)
Latitude: -43.01530 (43° 00' 55" S)
NAC:W80CK 7SWP3



More information about the postgis-users mailing list