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-
requ...@postgis.refractions.net wrote:
Message: 19
Date: Wed, 15 Feb 2012 10:04:01 +0100
From: Birgit Laggner <birgit.lagg...@vti.bund.de>
Subject: Re: [postgis-users] Landtype minus parcels ?
To: PostGIS Users Discussion <postgis-users@postgis.refractions.net>
Message-ID: <4f3b7501.2040...@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;
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
