Shapefile Reading: MultiPolygon being detected as (invalid) Polygon

[Robert Coup]

The Shapefile reader seems to detect multipolygons (2x exterior rings, no holes) as a polygon (1x ring, 1x hole) sometimes. Since the hole and ring don't overlap, the polygon is invalid.

In other cases it works fine, so I guess it's specific to some geometries.

The geometry was correctly detected in 1.6.1 but fails in 1.7.3, 1.8.0 and trunk, so I guess it's a regression.

Steps to reproduce:

• save test.csv and test.vrt
• run ogr2ogr -f "ESRI Shapefile" test.shp test.vrt
• see that under 1.6.1 you get a multipolygon, but under 1.7.3+ you get an invalid polygon geometry.

Result under 1.6.1:

$ ogrinfo -al test.shp 
INFO: Open of `test.shp'
      using driver `ESRI Shapefile' successful.

Layer name: test
Geometry: Polygon
Feature Count: 1
Extent: (175.524710, -40.172865) - (175.531271, -40.172035)
Layer SRS WKT:
GEOGCS["GCS_NZGD_2000",
    DATUM["New_Zealand_Geodetic_Datum_2000",
        SPHEROID["GRS_1980",6378137,298.257222101]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
id: String (80.0)
WKT: String (80.0)
OGRFeature(test):0
  id (String) = 1
  WKT (String) = MULTIPOLYGON(((175.5247097667 -40.17203475,175.5247128 -40.17207405,175.52473903
  MULTIPOLYGON (((175.524709766699999 -40.172034750000002,175.524757883299998 -4...

And under 1.7+ (including trunk at r22536):

$ ogrinfo test.shp -al
INFO: Open of `test.shp'
      using driver `ESRI Shapefile' successful.

Layer name: test
Geometry: Polygon
Feature Count: 1
Extent: (175.524710, -40.172865) - (175.531271, -40.172035)
Layer SRS WKT:
(unknown)
id: String (80.0)
WKT: String (80.0)
OGRFeature(test):0
  id (String) = 1
  WKT (String) = MULTIPOLYGON(((175.5247097667 -40.17203475,175.5247128 -40.17207405,175.52473903
  POLYGON ((175.524709766699999 -40.17203475,175.524757883299998 -40.17205056670...

[Robert Coup]

test data

[Robert Coup]

I couldn't find a simpler geometry that reproduced this, even playing with matching winding orders, etc.

The shapefiles produced by 1.6 and 1.8 are identical, so it is only the reader behaviour that's changed. I've attached the Shapefile now (test.zip).

[Robert Coup]

shapefile version of test.csv data

[warmerdam]

Even,

Is this an area you have done quite a bit of work in? Would you like to take this ticket?

[Even Rouault]

I've tried to compare quickly the related code between 1.6.1 and 1.7.3 and couldn't find a change that could explain a difference of behaviour, but whatever, let's move on and fix the bug...

So it is indeed a bug in OGRLinearRing::isClockwise(), which is revealed by the geometry part at the "south-east". The logic to select the points that will be used to compute the cross-product and thus deduce the orientation was too clever. In this particular situation trying to skip too close points to select others that are not too close from the origin result in selecting points whose orientation is the inverse of the expected orientation, hence the bug. As a fix, when the points are too close, I just fallback to computing the signed area with Green's formula, which works in your case and hopefully should work in most cases without introducing regressions in cases that worked (but I won't bet too much. It seems that each time I touch this code, someone comes with a new odd case that reveals a new bug...). At least, autotest still work and I've added your case for posterity.

Note that the bug is revealed by the optimization done in r15517, so you can workaround it by defining the configuration option/environmenet variable OGR_DEBUG_ORGANIZE_POLYGONS to YES (this will perform full "topological" analysis and disable the optimization).

For the record, this bug is linked to the changes done for #3356 and #3363.

Fixed in trunk (r22555) and in branches/1.8 (r22556)