[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