# [postgis-users] Getting dimensions of polygon sides

Obe, Regina robe.dnd at cityofboston.gov
Wed Jun 11 12:31:26 PDT 2008

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

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.

gives you the side from point 1 to point 2.  You would then need to loop

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