twpayne / go-geom

Package geom implements efficient geometry types for geospatial applications.
BSD 2-Clause "Simplified" License
847 stars 105 forks source link

The area of some 2D polygons is incorrect #224

Closed joosto closed 1 year ago

joosto commented 1 year ago

Hi there!

First of all, thank you for putting together this library, it's been a great help for the project I'm working on at the moment.

I think I've spotted an issue with area calculation of some 2D polygons with an interior ring. To caveat this, I'm entirely new to working with geospatial data, so let me know if you think otherwise. If you agree with the issue I'm happy to open a PR to fix this.

The problem

In my testing I found that the polygon I included in GeoJSON to the bottom of this issue, which has an exterior and an interior ring, has a larger area if you take into account the interior ring than if you consider the polygon without the cutout of the interior ring.

The issue is that the area of the interior ring is negative, so when subtracted from the area of the exterior ring the overall area is larger, which is incorrect.

I found on Wikipedia that the area of a polygon will be considered negative if it's negatively oriented. If that's the case, you can take the absolute value to get the area of the polygon.

If the polygon is negatively oriented, then the result A of the formulas is negative. In any case |A| is the sought area of the polygon.

Potential solution

Wrapping math.Abs around doubleArea on line 263 should fix the issue.

https://github.com/twpayne/go-geom/blob/3538962c3607662193a25804d1aef9e6c4f8eac3/flat.go#L258-L264

Example polygon that's affected by this issue

{"type":"Polygon","coordinates":[[[119918.67718197108,485786.741805515],[119918.69838706788,485785.29893093504],[119918.88923284179,485783.2831503406],[119919.24971930924,485781.3734634614],[119919.77984644996,485779.4425579076],[119920.43720410654,485777.6601834242],[119921.30661262499,485775.8778090019],[119922.2636510625,485773.72258849995],[119921.73365106489,485773.42658845807],[119923.16265109007,485770.86958843237],[119923.80765108918,485771.1265884876],[119926.17865111862,485768.61858840234],[119947.56065127277,485770.90258856124],[119948.17345583934,485771.039935658],[119949.33973553519,485771.35821677576],[119950.80288642865,485771.90990408685],[119952.11760171659,485772.5464664185],[119953.68677802975,485773.5861849084],[119955.19233908955,485774.7956532407],[119956.74031031624,485776.4082777353],[119957.60971880952,485777.63896488526],[119958.37310187577,485778.8908706995],[119959.11527985494,485780.4822764615],[119959.64540698321,485781.8614948073],[119959.98468834005,485783.049744464],[119960.21794427553,485784.15311904956],[119960.34517478317,485785.6808686479],[119960.34517478317,485787.7603053907],[119959.96348321895,485790.1368047347],[119959.41215098523,485791.7494291324],[119958.43671705364,485793.2347411709],[119957.80056448438,485794.041053401],[119956.82513055089,485795.1019906065],[119954.30172537726,485797.6694585515],[119952.90218973215,485798.8789268689],[119951.2693981506,485800.130832687],[119948.97924891827,485801.82833215175],[119937.8206509836,485810.0995887836],[119934.1136509639,485808.8275887822],[119933.91665095915,485809.31858880114],[119930.12665093935,485807.5255888168],[119930.46165094808,485807.0345888033],[119929.17369913645,485806.30548697076],[119926.96837026457,485804.6928623607],[119925.46280921216,485803.2712066501],[119924.2329142613,485801.95564451686],[119923.42712102513,485800.9583636061],[119922.62132778516,485799.81255143834],[119921.83673963061,485798.56064560346],[119920.84010063173,485796.45999000553],[119920.05551248007,485794.2108032379],[119919.46177010139,485792.4496474469],[119919.058873486,485790.5611792903],[119918.74079722205,485788.54539868433],[119918.67718197108,485786.741805515]],[[119938.367,485783.35],[119938.4305,485783.4683],[119938.455,485783.514],[119938.672,485784.816],[119938.71029786766,485784.8850877253],[119938.718,485784.952],[119938.789,485785.085],[119938.6211,485786.9466],[119938.598,485787.203],[119938.62227495036,485787.2467909786],[119938.604,485787.377],[119939.2457241494,485788.5346430877],[119939.665,485789.402],[119939.72748391585,485789.5147183282],[119939.732,485789.531],[119939.85304413277,485789.7493584703],[119940.054,485790.945],[119940.08762402345,485791.0056563091],[119940.0923,485791.0865],[119942.51977688179,485795.46556510986],[119942.67047688179,485795.5700651099],[119942.80947688178,485795.50606510986],[119942.80747688179,485795.47406510985],[119942.98337688179,485795.4380651099],[119943.67637688179,485795.2961651099],[119944.11147688178,485795.20706510986],[119944.22107688179,485795.2163651099],[119944.25247688178,485795.2190651099],[119944.27887688178,485795.2079651099],[119944.39247688178,485795.1600651099],[119944.55237688178,485795.16786510986],[119946.48530164218,485795.2620817077],[119946.50247688178,485795.2930651099],[119946.66847688178,485795.3090651099],[119946.79547688179,485795.2140651099],[119946.77847688178,485795.1760651099],[119947.74317688179,485794.7005651099],[119948.05657688179,485794.5461651099],[119948.66947688178,485794.2440651099],[119948.80147688178,485794.2190651099],[119948.93347688178,485794.1550651099],[119948.92347688178,485794.12006510986],[119949.90877688178,485793.9028651099],[119950.23447688179,485793.8310651099],[119950.40847688178,485793.84306510986],[119950.42047688179,485793.62606510986],[119950.37790957323,485793.54927549465],[119950.10547688178,485792.26706510986],[119950.13947688178,485792.1320651099],[119950.12147688179,485791.9930651099],[119950.09238926535,485791.9405922863],[119950.17087688179,485790.0332651099],[119950.17747688178,485789.8730651099],[119950.25047688179,485789.7220651099],[119949.38759609578,485788.16546485043],[119949.13847688178,485787.67406510987],[119949.118377559,485787.63780678593],[119949.14147688179,485787.5670651099],[119949.08047688179,485787.45206510986],[119949.05097688179,485787.33346510987],[119948.76157688178,485786.1695651099],[119948.75647688178,485786.1490651099],[119948.79147688179,485786.1440651099],[119948.78347688178,485785.9690651099],[119946.356,485781.59],[119946.236,485781.525],[119946.049,485781.587],[119945.5663,485781.2557],[119945.1974,485781.0026],[119945.165,485780.9803],[119945.013,485780.876],[119944.9559,485780.8985],[119944.861,485780.936],[119944.749,485780.769],[119943.288,485781.4139],[119943.0598,485781.5146],[119942.377,485781.816],[119942.216,485781.752],[119942.072,485781.841],[119942.08718697626,485781.8683966597],[119940.19,485782.794],[119940.192,485782.824],[119940.1716,485782.825],[119940.039,485782.832],[119939.916,485782.907],[119939.905,485782.935],[119939.8406,485782.9473],[119939.1416,485783.0809],[119938.634,485783.178],[119938.636,485783.208],[119938.449,485783.219],[119938.367,485783.35]]]}
twpayne commented 1 year ago

Thanks for your kind words!

The issue is a bit more subtle here. As you have identified, the orientation of the polygon determines whether its area is positive or negative. go-geom follows the convention that a polygon's outer ring has positive area and its inner rings have negative area, so the area can be calculated by a simple sum of all the ring areas: inner rings are subtracted without the need for a call to math.Abs.

In your case, the GeoJSON has incorrect polygon orientation for the inner rings. The GeoJSON spec states that holes are clockwise. So, the correct fix is either to fix the program that generated the GeoJSON or add a sanitization step when importing the data.

joosto commented 1 year ago

Ah that makes a lot of sense, thanks! That link to the GeoJSON spec is really useful.

The GeoJSON was generated by go-geom, but it's quite possible that the polygon I'm marshalling and calculating the area for is malformed. I'll have a look at that.

Thank you for your help 🙌