Simple conversion to Meshes.jl

[ErickChacon]

I want to convert ArchGDAL geometries to other types like Meshes.jl using GeoInterface.jl.

If I try to do it directly, it works without problems:

import ArchGDAL as AG
import GeoInterface as GI
using Meshes

mypoint = AG.createpoint(1.0, 2.0)
GI.coordinates(mypoint) |> Point

However it does not work if I try to do it using a custom function (it just hangs).

function topoint(geom)
    coords = GI.coordinates(geom)
    Point(coords)
end
topoint(mypoint)

Can anyone provide an insight of why it happens?

[rafaqz]

You really need an MWE showing what you ran exactly and the full stack trace of the error, or its very hard to know

(oops I missed that you said it just hangs... so there is no output at all? Try Ctrl-C and see what it was doing when you killed it. Also try just doing GI.coordinates(geom) in the function. If the problem is just in Point(coords) then this is certainly an issue for Meshes.jl, not ArchGDAL.jl, because coords is just vectors of numbers.).

But these days we usually use GeoInterface.convert(PackageName, object) to convert geometries between packages, rather than the old coordinates.

This works with GeometryBasics.jl, GeoInterface.jl, LibGEOS.jl, ArchGDAL.jl as sinks, and those plus Shapefile.jl, GeoJSON.jl, KML.jl and some others as sources (you can just write with those directly from any type so converting is redundant).

Meshes.jl has not implemented this interface as far as I know, and doesnt even depend on GeoInterface.jl. I think there is another package that implements some of it somehow for it though.

So this is not a problem with ArchGDAL.jl, but for Meshes.jl via GeoInterface.jl.

[ErickChacon]

My MWE is:

import ArchGDAL as AG
import GeoInterface as GI
using Meshes

mypoint = AG.createpoint(1.0, 5.0)

function topoint(geom)
    coords = GI.coordinates(geom)
    Meshes.Point(coords)
end

topoint(mypoint)

Unfortunately, I can not show the full stack trace of the error because julia hangs for ever. and Ctrl-C does not work, so I am forced to kill Julia. This same code work with other packages like Shapefile.jl and GeoJSON.jl, but does not work with ArchGDAL.

I am trying to check this basic function to help with the GeoInterface implementation for Meshes.jl which lives in GeoTables.jl (https://github.com/JuliaEarth/GeoTables.jl/blob/master/src/conversion.jl). I also open an issue in GeoTables (JuliaEarth/GeoTables.jl#13).

I want to use something like GeoInterface.convert, but inside this function I can see that other packages still call GI.coordinates. I tried to use GeoInterface.convert but Julia still hangs when I call Meshes.Point(GI.coordinates(geom)) inside any function, but only happens when geom is an ArchGDAL geometry.

I will continue checking what is the root of the problem, but happy to receive any suggestion.

[rafaqz]

Mostly they don't call coordinates? Im going through and changing that everywhere because its slow. (But yes ArchGDAL does, and that's not great)

You need to try it without the Meshes.Point - get the result of GI.coordinates(geom) from ArchGDAL.jl and inspect it, then work out why Meshes.jl cant handle it.

Running two functions together like that isnt possible to debug if you cant get a stack trace. Break it into smaller peices.

[ErickChacon]

I am closing this issue because it can be solved by making a constructor of Meshes.jl for Vec over AbstractVector type-stable (JuliaEarth/GeoTables.jl#13). Thank you so much for your help @rafaqz.