[postgis-users] Getting dimensions of polygon sides

[Obe, Regina]

 Bruce,

Thanks for the below.  I think it will take me a while to digest what
you are doing.  I think I learned a couple of things off the bat.

1) I may need the Right Hand Rule call since I was assuming things were
already ordered correctly and my simple spot check seems to suggest that
it is but you never know.

2) I'm confused between ST_Boundary and ST_ExteriorRing and when to use
one over the other.  I assume ST_Boundary takes into consideration holes
where as ST_ExteriorRing just gives you the outer ring.  Now I'm
thinking about it I probably should replace my ST_Boundary with
ST_ExteriorRing.

In case anyone is curious.  Attached is a snapshot of what I get when I
use my simple 2 point assumption.  It works in most cases but in others
where they used more than 2 points to describe each corner, I get extra
measures. 
Not the dimok is my ideal case and the dimalmost okay - see how I have
an extra measure.

 You think Simplifying before I extract would fix that.

My final query and updates look like this - seems to run fairly fast for
the 5000 parcel list I care about.

  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 ap.parcelid as pid, 2008 as pid_year,
ST_Boundary(ap.the_geom) As
the_boundary
			FROM dnd.rems_survey ap  ) p) As b
				CROSS JOIN generate_series(1, 1000) n
	WHERE n < b.nump;


UPDATE assessing.parcdimstime_2008 SET side_length =
CAST(ST_Length(the_geom) As numeric(8,2));



-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of
Bruce Rindahl
Sent: Wednesday, June 11, 2008 11:59 AM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Getting dimensions of polygon sides

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
>
>
>   

_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dimalmostok.png
Type: image/png
Size: 11422 bytes
Desc: dimalmostok.png
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20080611/1b318508/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dimok.png
Type: image/png
Size: 11662 bytes
Desc: dimok.png
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20080611/1b318508/attachment-0001.png>