[postgis-users] Landtype minus parcels ?

[Brian Hamlin]

Hi Birgit-

   thank you for the reply.. I admit that I have not compared your  
solution to my own carefully,
but instead have written a lot more code, and solved the problem.

   using psycopg2 and python's threading and Queue , I built a multi- 
core parallel version.
Since the logic loops through each large landtype geometry, and  
INSERTs the results
into one table, it is naturally parallelizable.  The psuedo-code goes  
something like this:

tSQL_primediff =
     INSERT into destination_table
       SELECT desired_fields,
          CASE when st_difference(
               landtype_poly, parcels
           ) is null then landtype_poly
          ELSE st_difference(
               landtype_poly, parcels )
          END as tgeom
FROM
     landtype_poly  INNER JOIN parcels
      ON st_intersects( landtype_poly, parcels)
WHERE
      landtype_poly id = tIND
GROUP BY
     landtype_poly id, landtype_poly

##---------------------------------------------
create a thread_UnitofWork subclass
     def run()
          get an id_start, id_end from queue msg
          SELECT landtype IDs >= start && <= end
          for each ID:
               cursor.execute(  tSQL_primediff )

          conn.commit()
          thread.done()

##----------------------------------------------------
      build a Queue
      SELECT landtype_id ORDER BY landtype_id
      count them and decide on a division per thread

      for the number of threads:
          build a queue entry
          fire the thread task

     wait for all threads to complete
###-------------------------------------------------------------
## last trick

tSQL_pre =
     SELECT lt_id FROM landtypes
       LEFT JOIN result_table  ON (lt_id = res.id )
     WHERE
        id_lt == NULL

tSQL_test =
    SELECT
       parcels_id
    FROM
       landtypes INNER JOIN parcels
         ON st_intersect( landtype, parcel)
    WHERE
         landtype_id = tIND

tSQL_alt =
     INSERT into result_table
      SELECT desired_fields, geom
      FROM landtypes
      WHERE landtypes_id = tIND


execute tSQL_pre, get a list of IDs
for each ID
    execute tSQL_test
    if  len(result_list) > 0: continue

    execute tSQL_alt  ## copy whole landtype into result

##----
the reason for this last step is that the way the program ran,
if a landtype was not at all touching any parcels, it was not copied
but you do not want to copy landtypes that were completely covered..

So you get the difference of landtype IDs in the input, and the ones  
that were
copied to the result, and check each one for the number of  
intersections of parcels.
Only if there were no intersections of parcels to that landtype, copy  
it to the result.

##-----
the final program uses 12 cores on a 16 core linux machine,  
processing sets of
landtypes in parallel.. It runs quite quickly on approximately  
850,000 parcels
against 14,000 landtype polys.

   thanks for the pointers and guidance.. perhaps I will get a chance  
to clean this up
at some later date..

==
Brian Hamlin
GeoCal
OSGeo California Chapter
415-717-4462 cell





On Feb 15, 2012, at 12:00 PM, postgis-users- 
request at postgis.refractions.net wrote:

>
> Message: 19
> Date: Wed, 15 Feb 2012 10:04:01 +0100
> From: Birgit Laggner <birgit.laggner at vti.bund.de>
> Subject: Re: [postgis-users] Landtype minus parcels ?
> To: PostGIS Users Discussion <postgis-users at postgis.refractions.net>
> Message-ID: <4F3B7501.2040000 at vti.bund.de>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> Hi Brian,
>
> you could try this query:
>
> SELECT st_difference(l.wkb_geometry,
>           case when st_union(p.wkb_geometry) is null
>             then st_collect(p.wkb_geometry)
>             else st_union(p.wkb_geometry)
>           end)
>         as tgeom
> from
>    big_landtypes l
>    inner join
>    small_parcels p
>    on st_intersects(l.wkb_geometry, p.wkb_geometry)
> where l.pkey = 1
> group by l.pkey, l.wkb_geometry;
>
> You are getting more than one result per landtype geometry because  
> they
> intersect with more then one parcel geometry. For each intersect case
> there is a result row. You are unioning all parcel geometries which  
> are
> anywhere intersecting any landtype geometry, which is probably  
> resulting
> in a big geometry. I would assume, the proposed query would do better,
> but you would have to try :-) I also added a CASE WHEN with st_collect
> as an alternative for st_union, because, sometimes st_union fails and
> results in a NULL geometry. In this case st_collect could take over  
> and
> as for the st_difference, st_union or st_collect have the same impact.
>
> Hope that helps,
>
> Birgit.
>
>
>
> Am 14.02.2012 20:42, schrieb Brian Hamlin:
>> SELECT distinct on (tgeom)
>>      st_difference(
>>        l.wkb_geometry,
>>
>>        (select st_union(p.wkb_geometry )
>>        from
>>          big_landtypes l,
>>          small_parcels p
>>        where st_intersects(
>>          l.wkb_geometry, p.wkb_geometry) AND
>>          l.pkey = 1)
>>        ) tgeom
>>      FROM
>>        big_landtypes l,
>>        small_parcels p
>>      WHERE
>>        l.pkey = 1
>>      ORDER BY
>>        tgeom;
>