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

Wed Apr 16 16:17:32 PDT 2008

Thanks Greg Troxel and Regina Obe for the input.

Based on Regina's ideas, I have something working.  Code is attached.  I 
have it working for 2 example parcels, one roughly square and convex, 
the other irregular and non-convex.  Map also attached.

I think I could modify this to take away some larger percentage of land 
on originally larger parcels, to apply a "net land use", to account for 
local roads and such, as Greg suggests.  I could also work from any 
corner or side (instead of the south-west corner), perhaps accounting 
for where the major roads are, as Greg suggests.   (We also have the 
option of telling model users that if there is not enough resolution in 
the model's output for a particular region, they could divide the 
polygon *before* our model runs, to allow other parts of our modelling 
system to take effect.)

I used a "view" for an interim table, instead of just an alias.  Is this 
going to slow things down?  The whole thing is quite slow on my two 
example parcels; I'm not sure how it's going to scale to 139000 parcels 
or 2.5 million parcels.  I'll add some more indices and then try it on 
the 139000 parcel set.

Finally, I don't know how to loop the last bit of code.  I could do it 
in java, python or a DOS Batch file, but if anyone can point me to how 
to do it in SQL, I would be much obliged.  Below are the commands that 
need to be executed in the first 2 runs through the final loop.  I can 
put the numbers "1" and "2" into a table easily enough, but for some 
reason I'm having a mental block about how to nest that into SQL.

UPDATE my_parcel_shards
    SET pseudoParcel = 1
    FROM bestPositionOption b
WHERE my_parcel_shards.parcel_no = b.parcel_no
    AND my_parcel_shards.shard_hpos <= b.shard_hpos
    AND my_parcel_shards.shard_vpos <= b.shard_vpos
    AND pseudoParcel IS NULL;

UPDATE my_parcel_shards
    SET pseudoParcel = 2
    FROM bestPositionOption b
WHERE my_parcel_shards.parcel_no = b.parcel_no
    AND my_parcel_shards.shard_hpos <= b.shard_hpos
    AND my_parcel_shards.shard_vpos <= b.shard_vpos
    AND pseudoParcel IS NULL;

Thanks so much,

--
John Abraham
jabraham at ucalgary.ca
-------------- next part --------------
-- create the pseduo_parcels, 20x20m.
drop table if exists my_parcel_shards;
Select p.parcel_no,
	p.hor_n as shard_hpos,
	p.ver_n as shard_vpos,
	0::int as initialPseudoParcel,
	null::int as pseudoParcel,
	ST_Intersection(p.the_geom,ST_Translate(
		p.boxref,p.hor_n*20,p.ver_n*20)) as the_shard
	into my_parcel_shards
	from (select parcel_no, hor.n as hor_n, ver.n as ver_n,
			the_geom,
			ST_SetSRID(ST_MakeBox2D(
				ST_MakePoint(ST_Xmin(the_geom),	ST_Ymin(the_geom)),
				ST_MakePoint(ST_Xmin(the_geom) + 20, ST_Ymin(the_geom) + 20)),
			26916) as BoxRef FROM pseudoparceltest
	cross join generate_series(0,99) As hor(n)
	cross join generate_series(0,99) as ver(n)) p
where
	ST_Intersects(p.the_geom, ST_Translate(p.boxref, p.hor_n*20, p.ver_n*20));

-- to restart
update my_parcel_shards set pseudoparcel=null;

-- add key field for QGIS
alter table my_parcel_shards add keyfield serial primary key;

--create view that tells us our options for taking bottom-left shards to make 2500to2900 sq m groups.
drop view positionoptions cascade;
create view positionoptions as 
SELECT f3.parcel_no, f3.shard_hpos, f3.shard_vpos, amt_used, 
		f3.shard_hpos*1.000001+f3.shard_vpos as compositePosition
		FROM (SELECT f1.parcel_no, f1.shard_hpos, f1.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_no = f2.parcel_no)
				WHERE   f1.shard_hpos >= f2.shard_hpos AND
					f1.shard_vpos >= f2.shard_vpos
					AND f1.pseudoParcel IS NULL
					AND f2.pseudoParcel IS NULL
				group by f1.parcel_no, f1.shard_hpos, f1.shard_vpos
				order by parcel_no, f1.shard_hpos, f1.shard_vpos
		) f3
			where f3.amt_used >= 6400 AND f3.amt_used <= 10000
			order by parcel_no, compositePosition;

-- pick the options that has the smallest "composite position"
--drop view bestPositionOption;
-- TODO this currently works its way up the diagonal.  It might be nice to 
-- count the vertical extent in npos and the horizontal extent in npos
-- and try to pick options that are more square, instead of options that are 
-- close to the diagonal.
create view bestPositionOption as select p.* from positionoptions p,
	(select min(compositeposition) as mincompositeposition, parcel_no from positionoptions group by parcel_no) f
where p.compositeposition = f.mincompositeposition and
	p.parcel_no=f.parcel_no;

UPDATE my_parcel_shards
    SET pseudoParcel = 1
    FROM bestPositionOption b
WHERE my_parcel_shards.parcel_no = b.parcel_no
    AND my_parcel_shards.shard_hpos <= b.shard_hpos
    AND my_parcel_shards.shard_vpos <= b.shard_vpos
    AND pseudoParcel IS NULL; 

UPDATE my_parcel_shards
    SET pseudoParcel = 2
    FROM bestPositionOption b
WHERE my_parcel_shards.parcel_no = b.parcel_no
    AND my_parcel_shards.shard_hpos <= b.shard_hpos
    AND my_parcel_shards.shard_vpos <= b.shard_vpos
    AND pseudoParcel IS NULL; 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: PseudoParcellingExample1.png
Type: image/png
Size: 9784 bytes
Desc: not available
Url : http://lists.refractions.net/pipermail/postgis-users/attachments/20080416/ff0e3b91/PseudoParcellingExample1.png