[postgis-users] Getting dimensions of polygon sides
Bruce Rindahl
rindahl at lrcwe.com
Wed Jun 11 08:58:41 PDT 2008
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
>
>
>
More information about the postgis-users
mailing list