On 2/11/2013 10:41 PM, Martin Feuchtwanger wrote:
> Do you mind posting a code snipet, Steve? Or maybe there are already
> good examples posted somewhere.
> I think yours is a classic case of GIS database analysis (not the exact
> example of soils-parcels-taxes, but the concept of deriving useful data
> based on the overlay of two independent layers).
> I'm new enough to Postgis to know only that it's going to involve
> ST_Intesection.
>
> Martin Feuchtwanger feumar at shaw.ca 604-254-0361
> http://members.shaw.ca/geomatics.developer
>
> On 11/02/2013 2:42 PM, Stephen Woodbridge wrote:
>>
>> I have done analysis where I had intersect a layer of parcel polygons
>> against a layer of soil types and then compute area of each soil types
>> contribution to the parcel and apply a tax rate for the soil to
>> compute the taxes on the parcel.
>
Hi Martin,
Conceptually the tax soils problem is not overly complicated. One of the
real challenges that I had was I did this 5 years ago on a much earlier
version of PostGIS and hard to deal with bugs and software that had a
lot less features. But that is not really important now.
I have not run any of the code below. So you might have to make changes
to get it to work, my goal is to present an algorithm and I'm using SQL
and PostGIS to present it.
We will store intermediate results in a new table. You could chain the
queries together and do it all at once, but if you are doing 1000's or
10s of 1000s of parcels on county sized region its better to save
intermediate results in case you need to come back to them as you can
see I make multiple queries. Five years ago when I did this job, it took
nearly a week to do this because of slower hardware and crashes and
having to rerun stuff multiple times to find and eliminate problematic
geometries. So depending on dataset sizes, run explain and make sure you
have the indexes you need, etc. and be patient.
So lets start out with two tables of polygons, parcels and soils. The
soils layer is a contiguous coverage of polygons, ie. it has no holes or
gaps and all the parcels are contained in the coverage. Our goal is to
chop the parcels into pieces based on the soils layer. We will store all
these pieces in a new table.
create table pieces (
gid serial not null primary key,
pid integer,
stype integer,
the_geom geometry -- Not using addgeometrycolumn because intersection
-- can return a collection
);
insert into pieces (pid, stype, the_geom)
select p.pid, -- parcel id
s.stype, -- soil type
st_intersection(p.the_geom, s.the_geom) as the_geom
from parcels p, soils s
where st_dwithin(p.the_geom, s.the_geom, 0.0);
At this point we should have all the intersection fragments in our new
table. Now we need to create a report of the areas of each soil type for
each polygon. If you had tax rate table for each soil type you could
also compute the taxes for each parcels.
-- compute the area by soil type and the tax by soil type by parcel
select p.pid,
p.stype,
sum(st_area(p.the_geom) as area),
sum(st_area(p.the_geom) * t.rate)
from pieces p, taxrates t
where p.stype=t.stype
group by p.pid, p.stype
order by p.pid, p.stype;
-- compute the total tax for each parcel
select p.pid, sum(st_area(p.the_geom) * t.rate) as tax
from pieces p, taxrates t
where p.stype=t.stype
group by p.pid
order by p.pid;
I also like to do some QA checks like do the some of the pieces of each
parcel equal the area of the original parcel? This is a sanity check and
they should be equal within numerical precision.
select a.pid, abs(st_area(a.the_geom) - b.area) as delta
from parcels a,
(select pid, sum(st_area) as area from pieces group by pid) as b
where a.pid=b.pid and delta > <limit>;
Some things to consider:
1. using an appropriate projection
2. if you have a lot a variance in parcel size you might want to change
the absolute limit check to a percentage check
3. I usually start any process like this be an isvalid() check of all
polygons and try to fix issues first.
Hope this helps,
-Steve