[encoding/geojson] RFC 7946 CRS WGS84

[TuSKan]

Hi,

Based on geojson RFC 7946, the coordinate system reference need to be WGS84. I think will be good to the geom ecosystem, save the CRS it are using, mainly when load from encoding packages. And also encoding/geojson always reproject the geom.T using your go-proj to WGS84 to compliment with RFC 7946. What do you think?

[twpayne]

No.

This would introduce a dependency go-proj, which requires cgo, and this is a pure Go library. There is no pure Go equivalent to PROJ, and, due to the complexity of PROJ, there is unlikely ever to be a pure Go version of PROJ. If you can use cgo, then consider using go-geos instead of go-geom.

[TuSKan]

I totally agree to not have dependency of go-proj. But what is your suggestion to have a geojson RFC 7946 complice ?

[twpayne]

But what is your suggestion to have a geojson RFC 7946 complice ?

What exactly do you mean by this?

Note that GeoJSON with non-EPSG:4326 projections is common and useful in the wild, and only the application using go-geom will know when it is handling such GeoJSON.

[TuSKan]

I understand. You are correct. In my application, I need a geojson RFC 7946 and I'm dealing with some difficulties to do that. I was able to do the reprojection to WGS84 using go-proj and geom.TransformInPlace. But can't figure out how to remove M from coordenates (change layout XYM to XY after reprojection).

[twpayne]

But can't figure out how to remove M from coordenates (change layout XYM to XY after reprojection).

A quick way to do it is:

pointXY := geom.NewPoint(geom.XY).SetCoords(pointXYM.Coords()))

You can do it more efficiently if you manipulate the flat coordinates directly, but this is only worthwhile if performance is a major concern.

[TuSKan]

Yes, performance is a major concern. After dealing with flat coordinate and stride, in a function like TransformInPlace, how to change g.layout ? There is no SetLayout and no SetFlatCoord

func ReProject(g geom.T, pj *proj.PJ, layout geom.Layout) {
    var (
        flatCoords = g.FlatCoords()
        stride  = g.Stride()
                strideIn = layout.Stride()
    )
    if strideIn > stride {
        panic("can't increase stride")
    } else if strideIn < stride {
        for i, n := 0, len(flatCoords); i < n; i += stride {
            c := flatCoords[i : i+stride]
            var d proj.Coord
            copy(d[:], c)
            o, err := pj.Forward(d)
            if err == nil {
                j := i - (stride-strideIn)*int(i/stride)
                copy(flatCoords[j:j+strideIn], o[:strideIn])
            }
        }
        // g.SetFlatCoords(flatCoords[:int(len(flatCoords)/stride)*strideIn])
        // g.SetLayout(layout)
    } else {
        for i, n := 0, len(flatCoords); i < n; i += stride {
            c := flatCoords[i : i+stride]
            var d proj.Coord
            copy(d[:], c)
            o, err := pj.Forward(d)
            if err == nil {
                copy(c, o[:])
            }
        }
    }
}

[twpayne]

There are multiple problems with the posted code:

• It combines re-projection with dropping extra dimensions.
• It does a lot of copying.
• It calls pj.Forward for every coordinate. This is a cgo call, which has a very high overhead.

If you don't need the extra dimensions, then don't request them from your data source. If you can't avoid requesting them, then create a new flat coords slice with the dimensions you want to keep and create a new geometry with those flat coords.

For the re-projection part, consider:

func ReProject(g geom.T, pj *proj.PJ) error {
    layout := g.Layout()
    return pj.TransFlatCoords(proj.DirectionFwd, g.FlatCoords(), g.Stride(), layout.ZIndex(), layout.MIndex())
}

This transforms all coordinates in-place with a single cgo call.

[TuSKan]

Thank you for the suggestion of ReProject. I didn't now about this function of flatCoords on go-proj.

This is my DropDimension func. Any suggestions are welcome.

func DropDimension(g geom.T, layout geom.Layout) geom.T {
    var (
        flatCoords = g.FlatCoords()
        ends       = g.Ends()
        endss      = g.Endss()
        stride     = g.Stride()
        strideNew  = layout.Stride()
    )

    if strideNew < stride {
        flatCoordsNew := make([]float64, int(len(flatCoords)/stride)*strideNew)
        for i := 0; i < len(flatCoords); i += stride {
            j := int(i/stride) * strideNew
            copy(flatCoordsNew[j:j+strideNew], flatCoords[i:i+strideNew])
        }
        endsNew := make([]int, len(ends))
        for i := 0; i < len(ends); i += 1 {
            endsNew[i] = ends[i] / stride * strideNew
        }
        endssNew := make([][]int, len(g.Endss()))
        for i := 0; i < len(endss); i += 1 {
            endssNew[i] = make([]int, len(endss[i]))
            for j := 0; j < len(endss[i]); j += 1 {
                endssNew[i][j] = endss[i][j] / stride * strideNew
            }
        }

        switch gg := g.(type) {
        case *geom.Point:
            return geom.NewPointFlat(layout, flatCoordsNew)
        case *geom.LineString:
            return geom.NewLineStringFlat(layout, flatCoordsNew)
        case *geom.LinearRing:
            return geom.NewLinearRingFlat(layout, flatCoordsNew)
        case *geom.Polygon:
            return geom.NewPolygonFlat(layout, flatCoordsNew, endsNew)
        case *geom.MultiPoint:
            return geom.NewMultiPointFlat(layout, flatCoordsNew)
        case *geom.MultiLineString:
            return geom.NewMultiLineStringFlat(layout, flatCoordsNew, endsNew)
        case *geom.MultiPolygon:
            return geom.NewMultiPolygonFlat(layout, flatCoordsNew, endssNew)
        case *geom.GeometryCollection:
            gColl := geom.NewGeometryCollection()
            for _, gi := range gg.Geoms() {
                gColl.MustPush(DropDimension(gi, layout))
            }
            return gColl
        }
    }
    return g
}   

[twpayne]

Looks good. You can avoid some multiplications with code like (warning: untested):

        for i, j := 0, 0; i < len(flatCoords); i, j = i+stride, j+strideNew {
            copy(flatCoordsNew[j:j+strideNew], flatCoords[i:i+strideNew])
        }

[TuSKan]

Thank you very much for your time helping me here. II think we can close this issue/helping for now

[TuSKan]

I'm sorry for using this place to make questions What you recommend for fast geometry search in go? (working with go-geom)

[twpayne]

What do you mean by "fast geometry search"?

[TuSKan]

Find a point inside a polygon, find the closest point with another

[twpayne]

Have you considered using Tile38 instead? Otherwise, you should be able to adapt github.com/tidwall/rtree.