[postgis-users] Cut off 1-acre of a parcel

[Obe, Regina]

John,

Very interesting problem.  First I am going to assume you have chosen a
Feet based spatial ref system and no parcel you have is greater than 
250 x 250 acres across. 

I'm going to break up your parcels into ~ 1/4 acres.  

I'm sure I'll screw up in my logic here, but hopefully it will be enough
to get you going.  

My idea is similar to yours except I am starting in south east corner
since that is the standard PostGIS origin

One way that comes to mind is to

1) First create a new table with your parcels clipped up into say  1/4th
acre units (I just need it to be lower than one acre to deal with
shards) with a parcel identifier denoting the parcel this slice is
formed from, the geometry, 2 fields -start_year (which you set to year
you started), ~year_used - which you update with your forecasting year. 

To do part 1 I would use a combo of ST_Translate, ST_Intersection,
generate_series something of the form.

SELECT p.parcel_id, hor.n As shard_hpos, ver.n As shard_vpos, 2008::int
As year_started, null::int As year_used,
ST_Intersection(p.the_geom, ST_translate(p.boxref, hor.n*width,
ver.n*height)) As the_shard
INTO my_parcel_shards
FROM (SELECT parcel_id,
ST_SetSRID(ST_MakeBox2D(ST_MakePoint(ST_Xmin(the_geom),
ST_Ymin(the_geom)),
	ST_MakePoint(ST_Xmin(the_geom) + 52, ST_Ymin(the_geom) + 52)), 
			your_srid_goes_here) As boxref 
		As boxref FROM parcels) p CROSS JOIN generate_series(0,
1000) As hor(n)
	CROSS JOIN generate_series(0,1000) As ver(n)
WHERE ST_Intersects(p.the_geom, ST_translate(p.boxref, hor.n*width,
ver.n*height))

2) Next continuously update the year_used field for each forecast year.
(This part I am sure I probably screwed up on, but hopefully
you get the idea).  The year_used marker will denote the year that shard
was used.

UPDATE my_parcel_shards
	SET year_used = whatever_forecast_year_we_are_up_to
	FROM (SELECT f3.parcel_id, MIN(shard_hpos) as hpos,
MIN(shard_vpos) As vpos
			FROM (SELECT parcel_id, shard_hpos, shard_vpos,
SUM(ST_Area(f2.the_shard)) As amt_used 
					FROM my_parcel_shards f1 INNER
JOIN my_parcel_shards f2 
						ON (f1.parcel_id =
f2.parcel_id)
						WHERE   f1.shard_hpos >=
f2.shard_hpos AND
	
f1.shard_vpos >= f2.shard_hpos
						AND f1.year_used IS NULL
AND f2.year_used IS NULL
					) f3   
				f3.amt_used >= 43560 AND f3.amt_used <=
46264
				GROUP BY f3.parcel_id) As forecast
WHERE my_parcel_shards.parcel_id = forecast.parcel_id
	AND my_parcel_shards.shard_hpos <= forecast.hpos
	AND my_parcel_shards.shard_vpos <= forecast.vpos
	AND year_used IS NULL;

Hope that helps,
Regina
 

-----Original Message-----
From: postgis-users-bounces at postgis.refractions.net
[mailto:postgis-users-bounces at postgis.refractions.net] On Behalf Of John
Abraham
Sent: Monday, April 14, 2008 8:11 PM
To: postgis-users at postgis.refractions.net
Subject: [postgis-users] Cut off 1-acre of a parcel

Hi there.  I'm a fairly new user of PostGIS, but I think I finally have 
it figured out and I'm finding it *very* useful and productive.

Now I'm ready to tackle one of my larger problems.  I'm working with a 
cadastral layer and forecasting future development patterns.  I have a 
need to simulate subdivision, where a large vacant parcel (say, 640 
acres) will realistically be developed in the future in smaller chunks 
(say, 1 acre).  Right now I ignore the spatial part of it, and end up 
having multiple future "pseudo-parcels" associated with one current-year

parcel.  So, for instance, if the simulation has two 1-acre development 
events in year 2010 and one 1-acre development event in year 2012, all 
occurring on the same original 640 acre parcel, then in year 2013 I have

4 records associated with the original 640 acre parcel: a 637 acre 
record that is still vacant, and three 1 acre records representing the 
development that has occurred.  All four of these records are still 
associated with the same polygon.

For the general predictive power of the model and for the mathematics 
behind it, this is ok, at least for now.  But for visualization, this 
isn't ok.  For visualization, I'd like a separate (non-overlapping) 
polygon for each of these records.  I need an algorithm that takes off 
an arbitrary 1acre chunk of the large parcel each time a development 
event is simulated for it.

The algorithm I'd like to use (unless someone suggests a better one) is.
1) find the east edge of the bounding box for the polygon.
2) draw a north-south line that is sqrt(43560) ft to the west of the 
east edge of the bounding box, and consider the areas on either side.
3) if the area on the east of the line is less than 1-acre, move the 
line west until the area on the east of the line is exactly one acre, 
done (stop).
4) else, draw an east-west line through the polygon far enough south so 
that the area north-east of the intersection of the two lines is exactly

1 acre.

I'm pretty sure this algorithm will grab 1-acre chunks of the original 
parcel starting in the north-east, and working down the east side, and 
when the east side is completely separated, it will start again at the 
new north-east corner and work down that side.   So the parcel will be 
slowly "subdivided" as development occurs from the east to the west in 
strips.  Since there are 43560 sqft in an acre, most of the chunks will 
be square.  (If the original polygon is non-convex, then the subdivided 
pieces may end up being multi-polygons, but that's ok, I can deal with 
multi-polygons, so long as the subdivision algorithm can handle it.  
There will also be some rather long-and-thin slivers, I think; it would 
actually be better to spiral around the parcel, working towards the 
inside in a circle, instead of always coming from the east side.  But...

I think this algorithm will work, and once it's working I can "tweak" it

later.)

Can anyone give me some pointers as to how to code this up in PostGIS?  
I have a few clues as to how to draw a line sqrt(43560) ft from the east

side, but I don't have any clues about how to draw a line "far enough 
south so that the resulting area is exactly 1 acre".  I could try some 
sort of trial-and-error bisection algorithm, I suppose, but I expect 
there might be a much easier way.

Any and all comments appreciated, thanks

--
John Abraham
jabraham at ucalgary.ca

_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
-----------------------------------------
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.