[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
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>
More information about the postgis-users
mailing list