I'm writing a python tool to update some river catchment boundaries and keep running into geometry errors, such as
WARNING: TopologyException: Input geom 1 is invalid: Self-intersection at or near point 1590337.0488225699 100863.34768028426 at 1590337.0488225699 100863.34768028426
If I call
.Buffer(0) on the geometry, the problem goes away, but then after writing out to shapefile and reading the feature back in, the problem returns.
I think I'm running into a floating point precision error - the coordinate system is in metres, so there is no need to have so many decimal places in the geometry coordinates. I found that if I read each point from the geometry, round to 2dp, then add the point to a new geometry, the problem goes away.
I don't believe there is a method in OGR or shapely to do this rounding (correct me if I'm wrong!), so I'm happy to do it manually. The issue is that not all of our data is in metres - sometimes the inputs are in WGS84, with map units in degrees. Rounding to 2dp is not going to help me here. What is the best way to decide how much to round my coordinates by?
I've tried using a percentage of
max(height, width), but it feels a bit clunky. Is it better to just have a lookup table of coordinate systems and round to, say 8dp if WGS84, or 2dp if anything else?