I am trying to use the Shapely and Fiona libraries to eliminate polygons from one layer (some triangles) from another layer (which represents the world's oceans).

I am looking for suggestions to improve the performance of how I've implemented this operation. I think I indexing the layers will help - using rtree or similar, but I haven't managed to get it to work.

Input Data:

  1. layer with polygon (7,3 MB) of the oceans worldwide
  2. layer with many triangles (300.000 aprox)

My goal is to find the fastest way to see whether or not the triangles are contained in the ocean polygons (1). For that, i have checked the prep method in shapely, but there are other methods such as str-packed-r-tree or using the library rtree in itself.

I was wondering what was the difference between them, how do they work, etc, to be able to guess what's the best to use in this case, and also to understand them better.

Below is the code of the function I wrote:

def clean_ocean(pathSourceOceans, pathSourceTriangles, pathOutput) :

    listGpkg = et1.fichiers_a_traiter(pathSourceTriangles,['gpkg'])
    print "list shapes", listGpkg

    for i in listGpkg :

        cheminOutput =  pathOutput+'netO_'+str(i)
        with fiona.open(pathSourceTriangles+i, 'r') as layerTriangles: 

            with fiona.open(pathSourceOceans, 'r') as layerOcean: 

                polygonOcean = asShape(layerOcean[1]['geometry'])
                prepPolygonOcean = prep(polygonOcean)

                # Create new schema 
                layerOutput_schema = {'geometry': 'Polygon','properties':{}}

                with fiona.open(
                        cheminOutput, 'w',
                        ) as layerOutput:

                    #do stuff:

                    for f in layerTriangles:      
                            # object.contains(other) "returns True if no points of other lie in the exterior of the object and at least one point of the interior of other lies in the interior of object.""
                            if not prepPolygonOcean.contains(asShape(f['geometry'])) :

                        except Exception, e:
                            logging.exception("Error processing feature %s:", f['id'])
  • Sorry, I just wanted to clarify a point about a "wrong" edit that has been made to my question, in the beginning: I have actually managed to get it work, what I am looking for is an explicit comparison about the different methods I've seen (prep(), rtree, or rtree embedded in shapely). Thanks anyway for the review! – bet_bit May 2 '18 at 13:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.