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 <[EMAIL PROTECTED]> 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: [EMAIL PROTECTED] 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@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: [EMAIL PROTECTED]
   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
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

Reply via email to