Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

The area of some 2D polygons is incorrect #224

Closed
joosto opened this issue Jan 2, 2023 · 2 comments
Closed

The area of some 2D polygons is incorrect #224

joosto opened this issue Jan 2, 2023 · 2 comments

Comments

@joosto
Copy link

joosto commented Jan 2, 2023

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.

go-geom/flat.go

Lines 258 to 264 in 3538962

func doubleArea1(flatCoords []float64, offset, end, stride int) float64 {
var doubleArea float64
for i := offset + stride; i < end; i += stride {
doubleArea += (flatCoords[i+1] - flatCoords[i+1-stride]) * (flatCoords[i] + flatCoords[i-stride])
}
return doubleArea
}

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
Copy link
Owner

twpayne commented Jan 2, 2023

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
Copy link
Author

joosto commented Jan 3, 2023

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 🙌

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants