4

From https://s.geo.admin.ch/77395652c7 I got the groundwater resources data for Switzerland. (deeplink to ZIP with Shapefiles, 147.3 MB).

Some areas therein are self-touching*, which makes gdalwarp unhappy when I try to use them to mask a raster:

gdalwarp cutline GK500/data/GK500_V1_3_Vector/GK500_V1_3_DE/Shapes/PY_Basis_Flaechen.shp -cwhere 'H2_ID = 6 OR H2_ID = 11' -crop_to_cutline 'MeteoSchweiz/Globalstrahlung/msg_SIS_Y_ch02_lonlat_2016_2016/msg.SIS.Y_ch02.lonlat_20160101000000_20160101000000.nc' out.tiff

results in

Warning 1: Ring Self-intersection at or near point 628843.78999999911 271660.44999999925 0
ERROR 1: Cutline polygon is invalid.

From PostGIS, I remembered the trick/hack to make such almost-valid polygons valid by "buffering" them with an infinitely small buffer (i.e. a buffer of width 0). This seems to work with OGR's python bindings, too:

from osgeo import ogr

gk500 = ogr.Open(
    'GK500/data/GK500_V1_3_Vector/GK500_V1_3_DE/Shapes/PY_Basis_Flaechen.shp')
gk_lyr = gk500.GetLayer()

for feature in gk_lyr:
    geom = feature.GetGeometryRef()
    if not geom.IsValid():
        feature.SetGeometryDirectly(geom.Buffer(0))
        assert feature.GetGeometryRef().IsValid()  # Doesn't fail

Though, when I then run

gk_lyr.ResetReading()

assert all(feature.GetGeometryRef().IsValid() for feature in gk_lyr)  # fails

that subsequent assertion fails. Why is that?

When iterating over the layer anew, are the features read from disk again?

Or is ogr.Feature.SetGeometryDirectly() not quite doing what I think it should do?

Please note that this question is not about how to solve this problem. If you want to provide your solution in your answer, fine, but I mostly want to know what is going on, so that I can form the right conceptual mental model about OGR.


*As far as I can tell, none are self-overlapping.

  • I suspect you may not be saving the changes, does this answer help? – Steven Kay Nov 30 '17 at 18:25
  • Actually, I don't want to save the changes to the disk, if I can avoid it. I just want to modify the in-memory representation and then work further with that within the same script, if possible. (Probably using the Python function gdal.Warp() instead of the gdalwarp CLI command.) – das-g Nov 30 '17 at 21:05
6

There are a couple small reasons that, in combination, are causing your assertion to fail.

Firstly, in the OGR/GDAL python bindings, any changes to files are only saved when the variable is deconstructed by either 1) the python garbage collector or 2) setting the python object equal to None to reset the variable.

If you are calling your two code blocks in sequence (in the same python code) without closing the variable (or allowing it to be closed automatically by the python garbage collector at the end of the script), the second assertion will fail because your changes to the shapefile have not been saved to the layer.

Method

I've developed my own code that I use quite frequently to clean shapefiles before I perform operations on them. It checks for both sliver polygons (by checking the area/perimeter ratio) and self-intersection polygons and generates a new shapefile (saved in memory to a python object) only containing valid geometries.

from osgeo import gdal, osr, ogr

def removeSliver(ds, EPSG):
    lyr = ds.GetLayer()
    out_ds = ogr.GetDriverByName('Memory').CreateDataSource('')
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(EPSG)
    out_lyr = out_ds.CreateLayer('', srs)
    for feature in lyr:
        geom = feature.GetGeometryRef()
        perimeter = feature.GetGeometryRef().Boundary().Length()
        area = feature.GetGeometryRef().GetArea()
        sliver = float(perimeter/area)
        if sliver <1: 
            if geom.IsValid() == True:
                wkt = geom.ExportToWkt()
                new_poly = ogr.CreateGeometryFromWkt(wkt)
                new_feat = ogr.Feature(out_lyr.GetLayerDefn())
                new_feat.SetGeometry(new_poly)
                out_lyr.CreateFeature(new_feat)
                new_feat = None
            else:
                new_feat = ogr.Feature(out_lyr.GetLayerDefn())
                new_feat.SetGeometry(geom.Buffer(0))
                out_lyr.CreateFeature(new_feat)
                new_feat = None
    out_lyr = None
    return out_ds

With this function it is very easy to clean a shapefile, and then do something to that shapefile:

ds = ogr.Open('/path/to/file.shp')
ds_clean = removeSliver(ds, 4326)
gdal.Warp('out.tif', 'input.tif', cutlineDSName=ds_clean, cropToCutline=True)

An added benefit is that the cleaned shapefile is only ever stored in memory while the original shapefile is never altered which maintains the integrity of the original shapefile.

Some Additional Notes:

Self-intersecting polygons and other invalid geometries don't actually causes errors in the output of gdalwarp, but instead are simply skipped when the utility attempts to perform a clip involving the invalid geometry. Depending on your use-case, you may or may not want to fix the geometry before clipping. Buffering a polygon by 0 will fix the self-intersection error but does not always produce intended results.

  • 1
    You miss a geom = feature.GetGeometryRef() somewhere in the for loop of removeSliver function. – Antonio Falciano Nov 30 '17 at 20:44
  • So is my assumption that re-iterating over the features of the layer re-reads the shapefile from disk (and thus discards my modifications that I didn't intend to write to disk but to use later in the script) correct? – das-g Nov 30 '17 at 21:19
  • 1
    Yes you are correct (but not for the reason you might think). GDAL/OGR is built on C+ so its all about whats going on with the pointers behind the scenes. In your second code block where you re-iterate through the features in gk_lyr with a for loop, the features are being pulled from the variable gk_lyr which points to the variable gk500 (your dataset object). However, the pointer has not been updated, so that pointer is still pointing to the original version of the shapefile (which contains invalid geometries), causing your assertion to fail. – onakua Dec 1 '17 at 0:10
  • It doesn't seem like cutlineDSName can take an existing ogr.DataSource instead of a filename. (Trying so results in ERROR 1: Cannot open <osgeo.ogr.DataSource; proxy of <Swig Object of type 'OGRDataSourceShadow *' at 0x7f11605b6690> >.) So I guess I have to write the cleaned shapes to a file, anyway. Nonetheless, thanks for bringing up ogr.GetDriverByName('Memory')! Didn't know about that, yet and may be handy in other situations. – das-g Dec 1 '17 at 11:03
  • You write "Self-intersecting polygons and other invalid geometries don't actually causes errors in the output of gdalwarp, but instead are simply skipped when the utility attempts to perform a clip involving the invalid geometry." This doesn't seem to be true in GDAL 2.2.1: Both, the CLI command and the Python function are aborting with ERROR 1: Cutline polygon is invalid. without producing the output file. – das-g Dec 1 '17 at 11:09
5

This is the debugged version of your code. You were trying to make valid the geometries without updating both the features and the datasource. See the inline comments to understand where the issues were.

from osgeo import ogr

gk500 = ogr.Open(
    'GK500/data/GK500_V1_3_Vector/GK500_V1_3_DE/Shapes/PY_Basis_Flaechen.shp', 1) # <====== 1 = update mode
gk_lyr = gk500.GetLayer()

for feature in gk_lyr:
    geom = feature.GetGeometryRef()
    if not geom.IsValid():
        feature.SetGeometry(geom.Buffer(0)) # <====== SetGeometry
        gk_lyr.SetFeature(feature) # <====== SetFeature
        assert feature.GetGeometryRef().IsValid() # Doesn't fail

gk_lyr.ResetReading()
assert all(feature.GetGeometryRef().IsValid() for feature in gk_lyr)  # Doesn't fail

Alternatively, I'd use a bit of SpatiaLite SQL to calculate the buffer of the invalid features, simply selecting the zero buffered not valid features together with the valid ones, without reinventing the wheel from scratch. It creates a new valid shapefile.

from osgeo import gdal

ds = gdal.OpenEx("GK500/data/GK500_V1_3_Vector/GK500_V1_3_DE/Shapes/PY_Basis_Flaechen.shp")

sql = "select ST_Buffer(geometry, 0) as geometry, rowid as fid from {} where ST_IsValid(geometry) = 0 \
      UNION select geometry, rowid as fid from {} where ST_IsValid(geometry) = 1" \
      .format(ds.GetLayer().GetName(), ds.GetLayer().GetName())
ds2 = gdal.VectorTranslate("PY_Basis_Flaechen_valid.shp", ds, SQLStatement = sql, SQLDialect = "SQLite")

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.