6

I'm trying to unpack shapefiles in Python (v3.7.4) using the GDAL-OGR library (v2.3.3). Each shapefile has only one feature containing multiple polygons with holes. They cannot be produced differently.

First I extracted the geometry and area. By comparing the areas calculated by GDAL and by ArcGIS, for some shapefiles the areas are different (GDAL's is bigger). By visual comparison, it is clear that GDAL failed to identify some holes:

Geometries and areas different in ArcGIS and GDAL

So, I unpacked each piece (as OGR multipolygon with holes model) and reconstructed the geometry (using polygons, outer and inner rings as discussed in Get the outer ring inner ring for a given polygon using ogr python binding). No luck. From the get-go, the number of geometries identified within the feature was 86, instead of 76. Here is the code:

from pathlib import Path
from osgeo import ogr

my_dir = Path(r'D:\Python')
shp = r'test_shape_area.shp'


Driver = ogr.GetDriverByName('ESRI Shapefile')
DataSource = Driver.Open((my_dir / shp).as_posix(), 0)
Layer = DataSource.GetLayer()
print('Number of features: {}'.format(Layer.GetFeatureCount()))

Geom = ogr.Geometry(ogr.wkbGeometryCollection)


for feature in Layer:
    multipolygon = feature.GetGeometryRef()
    print('Number of polygons: {}'.format(multipolygon.GetGeometryCount()))

    for j in range(multipolygon.GetGeometryCount()):
        polygon = multipolygon.GetGeometryRef(j)
        if polygon.GetGeometryCount() > 1:
            print('Number of rings (when >1): {}'.format(
                polygon.GetGeometryCount()))

        newPolygon = ogr.Geometry(ogr.wkbPolygon)

        for i in range(polygon.GetGeometryCount()):
            ring = polygon.GetGeometryRef(i)
            newPolygon.AddGeometry(ring)

        Geom.AddGeometry(newPolygon)


Geom.CloseRings()

print('Number of polygons in the final geometry: {}'.format(
    Geom.GetGeometryCount()))
print('Area of the final geometry: {}'.format(Geom.GetArea()))

Output:

In [1]: runfile('D:/Python/shp_with_holes_reader.py', wdir='D:/Python')
Number of features: 1
Number of polygons: 86
Number of rings (when >1): 5
Number of rings (when >1): 11
Number of rings (when >1): 4
Number of polygons in final geometry: 86
Area of the final geometry: 792794.0

It works if I transform the shapefile from multi to single features, but the idea is not to use ArcGIS. I tried lptorres' code (https://gist.github.com/lptorres/7781024) but ended up in the same place (it is basically the same unpacking).

So, I managed with shapely (based on GEOS):

from pathlib import Path
import shapefile
from shapely.geometry import shape

my_dir = Path(r'D:\Python')
shp = r'test_shape_area.shp'

my_shape = shapefile.Reader((my_dir / shp).as_posix())
print('Number of features: {}'.format(len(my_shape.shapeRecords())))
feature = my_shape.shapeRecords()[0]

first = feature.shape.__geo_interface__
shp_geom = shape(first)

print('Number of polygons in the final geometry: {}'.format(len(shp_geom)))
print('Area of the final geometry: {}'.format(shp_geom.area))

Output:

In [2]: runfile('D:/Python/shp_with_holes_reader_WORKING.py', wdir='D:/Python')
Number of features: 1
Number of polygons in the final geometry: 76
Area of the final geometry: 792314.0

Here is the link to the codes and shapefile: https://github.com/thainehassumpcao/shp_with_holes_reader

The question is: how to make it work with GDAL/OGR? It should work, right?


Update

Based on @J.Monticolo's comment, I've checked the multipolygon validity and tried to automatically correct it by buffering by 0 (code extract):

for feature in Layer:
    multipolygon = feature.GetGeometryRef()
    print('Number of polygons: {}'.format(multipolygon.GetGeometryCount()))

    if multipolygon.IsValid() == False:
        multipolygon = multipolygon.Buffer(0)

No correct result (as said, the holes were removed). Output:

runfile('D:/Python/shp_with_holes_reader.py', wdir='D:/Python')
Number of features: 1
Number of polygons: 86
Number of rings (when >1): 2
Number of rings (when >1): 2
Number of rings (when >1): 5
Number of rings (when >1): 5
Number of rings (when >1): 15
Number of rings (when >1): 3
Number of polygons in the final geometry: 75
Area of the final geometry: 792554.0
  • 1
    They are non-valid geometries, by trying to correct automatically the polygons, the identified holes disappear too. – J. Monticolo Jan 23 at 14:25
  • What software creates those polygons? – user30184 Jan 23 at 15:07
  • @J.Monticolo I see, I got that by buffering by 0 (updated the answer). Is there any other automatic correction that could be tested? In SQL I saw ST_MakeValid(), but couldn't find it in Python. – Thaine Herman Assumpção Jan 24 at 15:55
  • @user30184 HEC-RAS link, a flow simulation software. Unfortunately, I have no ways of altering the shapefile generation procedure. – Thaine Herman Assumpção Jan 24 at 16:00

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.