Can't read crs from simple flatgeobuf written with GDF

[alex-s-gardner]

import GeoFormatTypes as GFT
import GeoInterface as GI
using DataFrames
using GeoDataFrames

outfile = joinpath(tempdir(), "test.fgb");
pt = GI.Point(1,1, crs=GFT.EPSG(4326))

df = DataFrame(geometry=pt, test=1.)
GeoDataFrames.write(outfile, df)
test = GeoDataFrames.read(outfile)

produces this error

ERROR: Failed to convert this SRS into WKT format (Failure.)
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] macro expansion
    @ ~/.julia/packages/ArchGDAL/aXTAq/src/utils.jl:210 [inlined]
  [3] toWKT(spref::ArchGDAL.ISpatialRef)
    @ ArchGDAL ~/.julia/packages/ArchGDAL/aXTAq/src/spatialref.jl:591
  [4] read(ds::ArchGDAL.Dataset, layer::Int64)
    @ GeoDataFrames ~/.julia/packages/GeoDataFrames/eP9Sg/src/io.jl:96
  [5] (::GeoDataFrames.var"#10#11")(ds::ArchGDAL.Dataset)
    @ GeoDataFrames ~/.julia/packages/GeoDataFrames/eP9Sg/src/io.jl:62
  [6] read(f::GeoDataFrames.var"#10#11", args::String; kwargs::@Kwargs{})
    @ ArchGDAL ~/.julia/packages/ArchGDAL/aXTAq/src/context.jl:268
  [7] read
    @ ~/.julia/packages/ArchGDAL/aXTAq/src/context.jl:265 [inlined]
  [8] read(fn::String; kwargs::@Kwargs{})
    @ GeoDataFrames ~/.julia/packages/GeoDataFrames/eP9Sg/src/io.jl:58
  [9] read(fn::String)
    @ GeoDataFrames ~/.julia/packages/GeoDataFrames/eP9Sg/src/io.jl:52
 [10] top-level scope
    @ ~/Documents/GitHub/Altim.jl/src/import GeoFormatTypes as GFT.jl:70

other file types DO NOT have this issue:

import GeoFormatTypes as GFT
import GeoInterface as GI
using DataFrames
using GeoDataFrames

pt = GI.Point(1,1, crs=GFT.EPSG(4326))
df = DataFrame(geometry=pt, test=1.)

fn = ["test.geojson", "test.fgb", "test.shp", "test.gpkg"]
for f in fn
    outfile = joinpath(tempdir(), f);
    GeoDataFrames.write(outfile, df)
    try
        test = GeoDataFrames.read(outfile)
    catch e
        println("causes read error: $f")
    end
end

[evetion]

Thanks for making an issue! I will look into it.

[evetion]

Can you post the versions used? I can't reproduce it. To be fair, there's the bug (maybe related) that the crs of a geometry is not set on the DataFrame.

This is what I get from metadata(test)

  "crs"                          => WellKnownText{CRS}(CRS(), "LOCAL_CS[\"Undefined SRS\",LOCAL_DATUM[\"unknown\",32767],UNIT[\"unknown\",0],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")`

You can set the crs with

GeoDataFrames.write(outfile, df; crs=GFT.EPSG(4326))

[alex-s-gardner]

Can you post the versions used?

[62cb38b5] GeoDataFrames v0.3.10
[c9ce4bd3] ArchGDAL v0.10.4 `https://github.com/asinghvi17/ArchGDAL.jl#patch-3

[alex-s-gardner]

@evetion a couple of clarifying questions:

1. Do I need to specify the crs of the DataFrame even when the geometries contain crs info? If so is there an easy way to pull this from the geometries by default?

2. Why does the OP work for geojson, .shp, .gpkg but not .fgb?

[evetion]

1. At the moment yes. I will make a PR to fix it (but will have to error/warn when geomtries have different crs'es.
2. I think there's a GDAL .fgb bug in that an empty spatial reference is written as an empty string, which can't be parsed back. This is handled correctly in other drivers.

I didn't reproduce at first, because in the unreleased master, we set a default local CRS instead of nothing (v0.3.10...master#diff-db104fa7ac6c3097a327fc13a1076a2f1dd30a49057cae4caebe58563a97030eR175). Also, you don't need asinghvi17s patched fork anymore, ArchGDAL 0.10.5 is out.

[alex-s-gardner]

@evetion, thanks for all your help with this.