0

I am using fiona and shapely together to buffer points and then merge them and write them to a shapefile. For almost all of my datasets, this is working fine. However, I am currently getting a MemoryError on the mapping(merged) call.

Here's the error that I'm getting:

Traceback (most recent call last):

File "C:\Working\Projects\USACE_Data\Toolbox\V4.1\Tools\ImportEHydro.py", line 149, in execute helper.ExcludeNewerDataWithinArea(lyrUSACE, fileTemp, fcSurveyJob, long(sordat.strftime("%Y%m%d")))

File "C:\Working\Projects\USACE_Data\Toolbox\V4.1\Tools\HelperClass.py", line 92, in ExcludeNewerDataWithinArea self.USACEMinimumBoundingPolygon(usaceLayer, newerDataPoly , "SORDAT")

File "C:\Working\Projects\USACE_Data\Toolbox\V4.1\Tools\HelperClass.py", line 609, in USACEMinimumBoundingPolygon firstBufferFC = self.BufferFeatureClass(tempCopy, "new_buffer", mean + (4*std), mergeType, mergeField)

File "C:\Working\Projects\USACE_Data\Toolbox\V4.1\Tools\HelperClass.py", line 540, in BufferFeatureClass 'geometry': mapping(merged)

File "C:\python\ArcGIS10.5\lib\site-packages\fiona\collection.py", line 353, in write self.writerecords([record])

File "C:\python\ArcGIS10.5\lib\site-packages\fiona\collection.py", line 347, in writerecords self.session.writerecs(records, self)

File "fiona\ogrext.pyx", line 1215, in fiona.ogrext.WritingSession.writerecs MemoryError

The relevant bits of code are below:

with collection(os.path.join(newShpDir, outputShp + ".shp"), "w", "ESRI Shapefile", schema) as output:
     # Merge all of the features into one feature
     for i in range(len(features)):
         try:
             merged = cascaded_union(features[i])
         except:
             arcpy.AddWarning("  Could not merge features for " + (str(featureList[i]) if len(featureList > 0) else "(unknown)") + "...skipping...")
             continue

         # write the output to a shapefile
         dictProperties = {}
         if mergeAll:
             dictProperties = {mergeField : str(bufferDist)}
         else:
             dictProperties = {mergeField : featureList[i]}

         output.write({
             'properties': dictProperties,
             'geometry': mapping(merged)
         })

The error is happening at the 'geometry':mapping(merged) It's trying to merge about 350k shapes into one. Anyone know if there is a limit to how many it can handle or if this is normal? I can handle it with an exception and then union some, merge, union some, merge, etc...but I'd rather not if I don't have to...

  • I was initially assuming that it was a shapely problem...I missed that it was in fiona that it was getting a memoryerror and not the shapely library. It, theoretically, should just be writing a single feature geometry, but that one geometry is created from a buffer of 350k shapes that means there could be a LOT of vertices. I'm not really sure what to do because I can't split them up as I need them split up geographically... – William Winner Oct 4 at 18:11
1

Welcome to the site William. Yes, there are some geometry limitations that govern how much information can be stored in .shp format. It is likely that your 350k features are going over the 2GB file size limit.

See this page for more details on the constraints of the ESRI shapefile format. https://desktop.arcgis.com/en/arcmap/latest/manage-data/shapefiles/geoprocessing-considerations-for-shapefile-output.htm

  • 2
    I did some checking of the output from the 350k features and it was creating a shapely MultiPolygon with around 9m vertices. According to the ESRI documentation, the 2GB limit should allow around 70 million points so I'm not sure it was reaching that limit. But, I did add a try except checking for a MemoryError and, if it happened, I simplify the geometry until I don't get a Memory Error. That allowed my routine to run. – William Winner Oct 11 at 12:55
0

Just in case anyone wanted to see what I did, here are the changes to the code that allowed it to work:

with collection(os.path.join(newShpDir, outputShp + ".shp"), "w", "ESRI Shapefile", schema) as output:
    # Merge all of the features into one feature
    for i in range(len(features)):
        try:
            merged = cascaded_union(features[i])
        except:
            arcpy.AddWarning("  Could not merge features for " + (str(featureList[i]) if len(featureList > 0) else "(unknown)") + "...skipping...")
            continue

    # sub-routine to simplify the geometry
    def simplifyGeometry(geom, dist):
        newGeom = geom
        # if it's a MultiPolygon (shapely type) we will simplify the individual polygons
        if geom.geom_type == 'MultiPolygon':
            arcpy.AddMessage("Simplify type = MultiPolygon")

            # create a list of polygons
            geomList = list(geom)
            newList = []
            arcpy.AddMessage("Simplifying " + str(len(geomList)) + " features using tolderance of " + str(dist))

            # go through each polygon and simplify them.  We go through each to ensure we don't have any
            # bad geometries created
            for polygon in geomList:
                try:
                    newPoly = polygon.simplify(dist)

                    #fix holes
                    newPoly = newPoly.buffer(0.001, 1, join_style=JOIN_STYLE.mitre).buffer(-0.001, 1, join_style=JOIN_STYLE.mitre)
                except:
                    # got an error likely saying it couldn't handle a null geometry
                    #geomList.remove(polygon)
                    pass

                # make sure the new geometry is valid and add it to the new list
                if newPoly.is_valid:
                    newList.append(newPoly)

            # create another multipolygon
            newGeom = cascaded_union(newList)
        else:
            newGeom = geom.simplify(dist)
        return newGeom

    # sub-routine to write the output
    def writeOutput(geom, factor):
        try:
            # write the output to a shapefile
            dictProperties = {}
            if mergeAll:
                dictProperties = {mergeField : str(bufferDist)}
            else:
                dictProperties = {mergeField : featureList[i]}

            output.write({
                'properties': dictProperties,
                'geometry': geom
            })
        except MemoryError:
            arcpy.AddMessage("Memory Error writing to file; simplifying geometry")

            # simplify the geometry and increase the factor for next time
            newGeom = simplifyGeometry(shape(geom), bufferDist * factor)
            factor += 2

            # try writing the output again
            writeOutput(mapping(newGeom), factor)

        arcpy.AddMessage("Writing output to shapefile")

        writeOutput(geometry, 1)

Interestingly, I don't have to close and open the fiona collection because it never actually writes anything. I also would point out that collection has been deprecated but still works so if you are starting new, use fiona's open instead of collection.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.