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
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?
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.