[postgis-users] Sample usages of analytic functions

Stephen Woodbridge woodbri at swoodbridge.com
Tue Feb 12 21:30:37 PST 2013


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


More information about the postgis-users mailing list