Incorrect centroids on geometries with holes

[jerluc]

Observed

For geometries with holes, the planar.CentroidArea() function is returning incorrect centroids (the planar area returned is correct though).

To reproduce

Note that this bug has been observed on other geometries, but this is the one I've chosen as an example:

package main

import (
    "fmt"
    "github.com/paulmach/orb"
    "github.com/paulmach/orb/planar"
)

func main() {
    geom := orb.Polygon{
        {
            {
                -102.2690493,
                40.9939916,
            },
            {
                -102.268951,
                40.9940453,
            },
            {
                -102.2690464,
                40.9941449,
            },
            {
                -102.2689287,
                40.9942091,
            },
            {
                -102.2688078,
                40.9940829,
            },
            {
                -102.2684138,
                40.9942979,
            },
            {
                -102.2683195,
                40.9941995,
            },
            {
                -102.2684496,
                40.9941286,
            },
            {
                -102.2683239,
                40.9939973,
            },
            {
                -102.2682875,
                40.9940171,
            },
            {
                -102.2682667,
                40.9939954,
            },
            {
                -102.2680722,
                40.9941015,
            },
            {
                -102.2679264,
                40.9939493,
            },
            {
                -102.2680433,
                40.9938855,
            },
            {
                -102.268059,
                40.9939019,
            },
            {
                -102.2682353,
                40.9938057,
            },
            {
                -102.2682124,
                40.9937818,
            },
            {
                -102.2682594,
                40.9937561,
            },
            {
                -102.2682808,
                40.9937785,
            },
            {
                -102.2683036,
                40.993766,
            },
            {
                -102.2683272,
                40.9937907,
            },
            {
                -102.2685813,
                40.9936521,
            },
            {
                -102.2688361,
                40.9939181,
            },
            {
                -102.2689632,
                40.9938488,
            },
            {
                -102.2690246,
                40.993913,
            },
            {
                -102.2689913,
                40.9939311,
            },
            {
                -102https://github.com/paulmach/orb/blob/master/planar/area.go#L197.2690493,
                40.9939916,
            },
        },
        {
            {
                -102.2687189,
                40.9939914,
            },
            {
                -102.2685563,
                40.9938227,
            },
            {
                -102.2684289,
                40.9938926,
            },
            {
                -102.2684666,
                40.9939317,
            },
            {
                -102.2684333,
                40.99395,
            },
            {
                -102.2685582,
                40.9940796,
            },
            {
                -102.2687189,
                40.9939914,
            },
        },
    }
    centroid, _ := planar.CentroidArea(geom)
    fmt.Println(centroid)
}

Running this code produces the following output:

[-115.88899555567677 46.45368396555677]

Expected

Using other tools (in this case, shapely) produces the following, which can be visually verified as being correct/expected:

[-102.2685248962348 40.99397171799407]

Possible causes

One bit of code I ran into while trying to debug is in the implementation for planar.CentroidArea(). I have read in some other implementations (namely libgeos) that they use the "signed" area, rather than the absolute value, as orb uses here. I'm not 100% sure that this is the cause, but wanted to ask about this logic.

[jerluc]

I think I may have found the root of the issue, creating a PR now!

[paulmach]

This was fixed a month about via https://github.com/paulmach/orb/pull/16 and https://github.com/paulmach/orb/commit/f3ee61f36a55c286ab302585a71ecefaac5661f5