1

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?

jon_two
  • 1,854
  • 12
  • 19
  • https://gis.stackexchange.com/questions/8650/measuring-accuracy-of-latitude-and-longitude may help you – nmtoken Sep 16 '19 at 17:23
  • @nmtoken Not really, I just want to find out if the given geometry is in metres or degrees. – jon_two Sep 16 '19 at 19:28
  • You could start by checking to see if the coordinate system is a geographic coordinate system (GCS) or a projected coordinate system (PCS). GCS will be in degrees, PCS will not be in degrees. There is no guarantee that an arbitrary PCS will be in metres, but most of them are. You may be in the position of knowing that for you particular case all your PCS are in metres, so just distinguishing between geographic and projected coordinate system may be enough for you. ogrinfo would give you this information. – Son of a Beach Sep 17 '19 at 00:37
  • If the number in front of the decimal is greater than 180 it's unlikely to be in degrees, – nmtoken Sep 17 '19 at 06:05
  • 1
    Ok, so maybe just srs.IsProjected() would be enough in most cases? The majority of our data is only in four or five coordinate systems, but users do like to surprise us. – jon_two Sep 17 '19 at 06:39
  • 1
    I was also wondering if someone would come on to tell me rounding the coordinates was a ridiculous idea, but nobody has yet. Maybe I should just go with it. – jon_two Sep 17 '19 at 06:40

0 Answers0