[postgis-users] Getting dimensions of polygon sides

[Bruce Rindahl]

Regina
I don't know if this is exactly what you are looking for but it might 
give you ideas.
I have a polygon table of all the sections in Colorado (PLSS).  I needed 
to get the "North" side of the section as well as the "South, West and 
East".
First I made sure all the polygons were digitized in the same 
direction.  This was done via forceRHR(the_geom).
Next I added 4 attributes called ne,se,sw and nw.  These attributes hold 
the index number of the point closest to the corners of the section.  I 
then created a function called line_locate_vertex(line,point)  which 
returns the index needed above:
-----------------------------------------------------------------------------------------------------
-- Function: line_locate_vertex(geometry, geometry)
-- DROP FUNCTION line_locate_vertex(geometry, geometry);

CREATE OR REPLACE FUNCTION line_locate_vertex(geometry, geometry)
  RETURNS integer AS
$BODY$
  DECLARE    
    dist double precision;
    n_dist double precision;
    i integer;
    p integer;
  BEGIN
    dist = 9999999999999;
   IF (GeometryType($1) != 'LINESTRING') THEN
    raise notice 'First Geometry type must be a LINESTRING';
    RETURN -1;
   END IF;
   IF (GeometryType($2) != 'POINT') THEN
    raise notice 'Second Geometry type must be a POINT';
    RETURN -1;
   END IF;
    FOR i IN 1 .. numPoints($1)
    LOOP
        n_dist = distance($2,pointN($1,i));
        IF (n_dist < dist) THEN
            p = i;
            dist = n_dist;
        END IF;
    END LOOP;
    RETURN p;
  END;
$BODY$
  LANGUAGE 'plpgsql' VOLATILE;
ALTER FUNCTION line_locate_vertex(geometry, geometry) OWNER TO postgres;
----------------------------------------------------------------------------------------------------------------
Using this function, I can then update my polygon coverage with the 
index number of the section corners.
----------------------------------------------------------------------------------------------------------------
update sections set ne = 
line_locate_vertex(exteriorRing(geometryN(the_geom,1)),setSrid(MakePoint(xmax(the_geom),ymax(the_geom)),srid(the_geom)));
update sections set se = 
line_locate_vertex(exteriorRing(geometryN(the_geom,1)),setSrid(MakePoint(xmax(the_geom),ymin(the_geom)),srid(the_geom)));
update sections set sw = 
line_locate_vertex(exteriorRing(geometryN(the_geom,1)),setSrid(MakePoint(xmin(the_geom),ymin(the_geom)),srid(the_geom)));
update sections set nw = 
line_locate_vertex(exteriorRing(geometryN(the_geom,1)),setSrid(MakePoint(xmin(the_geom),ymax(the_geom)),srid(the_geom)));
----------------------------------------------------------------------------------------------------------------------------------
Now the "sides" of the polygons are all defined as the segment running 
from one of the corners to another.  For example the "North" section 
line might go from point 2 to point 5.  To return a particular side I 
created a function section_line(line,integer,integer) which returns a 
line segment from point $1 to point $2.
--------------------------------------------------------------------------------------------------------------------
-- Function: section_line(geometry, integer, integer)
-- DROP FUNCTION section_line(geometry, integer, integer);
CREATE OR REPLACE FUNCTION section_line(geometry, integer, integer)
  RETURNS geometry AS
$BODY$
  DECLARE    
    section_line GEOMETRY;
    i integer;
    last integer;  
  BEGIN
   IF (GeometryType($1) != 'LINESTRING') THEN
    raise notice 'Geometry type must be a LINESTRING';
    RETURN -1;
   END IF;
   SELECT into section_line $1;
    last = $3;
    IF ($2 > last) THEN
    FOR i in 2 .. ($3)
        LOOP
            SELECT into section_line 
AddPoint(section_line,PointN(section_line,i),-1);
            last = last + 1;
        END LOOP;
    ELSE
    FOR i in REVERSE NumPoints(section_line) .. ($3 + 1)
        LOOP
            SELECT into section_line RemovePoint(section_line,i - 1);
        END LOOP;
    END IF;
    FOR i in 0 .. ($2 - 2)
    LOOP
        SELECT into section_line RemovePoint(section_line,0);
    END LOOP;
    RETURN section_line;
  END;
$BODY$
  LANGUAGE 'plpgsql' VOLATILE;
ALTER FUNCTION section_line(geometry, integer, integer) OWNER TO postgres;
--------------------------------------------------------------------------------------------
Note this function creates a line segment of that side even if it has to 
union part of the start with part of the end.

All this now gets me to what I want - an arbitrary section line:

SELECT
section_line(exteriorRing(geometryN(the_geom,1)),nw,ne) as north,
section_line(exteriorRing(geometryN(the_geom,1)),se,sw) as south,
section_line(exteriorRing(geometryN(the_geom,1)),sw,nw) as west,
section_line(exteriorRing(geometryN(the_geom,1)),ne,se) as east
from sections WHERE .......

Using these functions I created a query where given an arbitrary point 
it returns the section, township, range and distance from the closest 
north/south line and the closest east/west line for a legal definition.

In your case section_line(exteriorRing(geometryN(the_geom,1)),1,2)  
gives you the side from point 1 to point 2.  You would then need to loop 
thorough all the points in your polygon to get your table.

Sorry about the long winded post - I hope this gives you ideas.  If you 
want I have all this in a single sql file.

Bruce Rindahl


Obe, Regina wrote:
> I'm trying to similuate ArcGIS annotations.  I guess ArcPad doesn't
> support ArcGIS annotations according to what I have been told but can
> support line strings and so forth.  
>
> Here is the the problem.  I have a set of parcel polygon geometries.
> I'm going to assume that each lines side is composed of 2 points and no
> polygon has more than 1000 sides.  I need to create a table that has a
> separate row for each side with the length  of that side as an attribute
> fields.
>
> I stupidly thought I could take the boundary and then figure out the
> length of the boundary forgetting that this just gives me the perimeter.
>
> So my second thought was that if I reconstitute points of the boundary
> grouping 2 at a time  - that would do the trick.
>
> My query is still running so haven't looked to see what the final result
> is.  I'm wondering if someone has done something similar and if they
> have an easier way.
>
> Below is the query I am testing right now.
>
>   INSERT INTO assessing.parcdimstime_2008(pid, pid_year, the_geom)
> 	SELECT b.pid, b.pid_year, ST_MakeLine(ST_PointN(the_boundary,n),
> ST_PointN(the_boundary, n + 1)) As the_side
> 	FROM (SELECT pid, pid_year, the_boundary,
> ST_NumPoints(the_boundary) As nump
> 		FROM (SELECT pid, pid_year, ST_Boundary(the_geom) As
> the_boundary
> 			FROM assessing.parceltime
> 			WHERE pid_year = 2008) p) As b
> 				CROSS JOIN generate_series(1, 1000) n
> 	WHERE n < b.nump;
>
> UPDATE assessing.parcdimstime_2008 SET sd_length = ST_Length(the_geom);
>
> Thanks,
> Regina
>
> -----------------------------------------
> 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.
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
>