Here is the code.
Am 22.01.2016 um 14:05 schrieb Kai Muehlbauer:
Am 22.01.2016 um 13:26 schrieb Even Rouault:
Le vendredi 22 janvier 2016 13:19:40, Kai Muehlbauer a écrit :
Am 22.01.2016 um 12:58 schrieb Ari Jolma:
22.01.2016, 13:49, Kai Muehlbauer kirjoitti:
This is very good news! If I understand correctly, I could use your
patch and self-compile gdal to have access to the improved
layer-intersection function via the python bindings. I'll give this a
try.
Attached is probably a bit better diff.
As I wrote, similar thing should be done for other layer methods.
Yes, you can use this with the trunk version.
Thanks Ari, I keep you posted. If I get this running, I'm also available
for testing the other changes.
Kai, just thinking: do you have spatial indices defined, especially on
the
method layer ? If not, that would be the first thing to do since for each
geometry in input layer, a spatial filter is set on method layer using
the
extent of the input geometry.
Thanks for the comment Even, yes, I have a spatial filter in place. But
with prepared geometries it takes only 50% of the time. Attached you
find my (stripped) python code:
Method 1 is more or less what I do. In my example this takes 35 seconds.
Method 2 I implemented with shapely prepared geometry. This runs only
~17 seconds on the same machine.
If I call layer.Intersection(method, dst, options) directly it takes ~55
seconds.
So with the prepared geometries we get 70% improvement. And 50%
improvement compared to the iteration over layer method.
And, without SpatialFilter this will run almost forever.
Thanks again for your comments, it's very much appreciated.
Cheers,
Kai
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
--
Kai Muehlbauer
Meteorological Institute University of Bonn
Auf dem Huegel 20 | +49 228 739083
D-53121 Bonn | [email protected]
# method 1
for feature in layer:
trg = feature.GetGeometryRef()
method.ResetReading()
method.SetSpatialFilter(trg)
for feature in method:
geom = feature.GetGeometryRef()
if not trg.Intersects(geom):
geom = trg.Intersection(geom)
new = ogr.Feature(defn)
new.SetGeometry(geom)
dst.CreateFeature(new)
# method 2
from shapely import wkb
from shapely import prepared
for feature in layer:
trg = feature.GetGeometryRef()
# get shapely objects
trg_shape = wkb.loads(trg.ExportToWkb())
trg_prep = prepared.prep(trg_shape)
method.ResetReading()
method.SetSpatialFilter(trg)
for feature in method:
geom = feature.GetGeometryRef()
# if not trg.Intersects(geom):
# prepared contains test
if not trg_prep.contains(wkb.loads(geom.ExportToWkb())):
geom = trg.Intersection(geom)
new = ogr.Feature(defn)
new.SetGeometry(geom)
dst.CreateFeature(new)
_______________________________________________
gdal-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
