yeesian / ArchGDAL.jl

A high level API for GDAL - Geospatial Data Abstraction Library
https://yeesian.github.io/ArchGDAL.jl/stable/
Other
142 stars 26 forks source link

Access geometry of layer converted to dataframe #265

Closed pwidas closed 3 years ago

pwidas commented 3 years ago

Thanks for the great work on developing this package. I'd like to access the numerical values of the geometry in a dataframe that is created with ArchGDAL and Tables.jl. With the code below (and data attached) I get the error "Incompatible geometry for operation". What's my error?

julia> import ArchGDAL as AG julia> import DataFrames as DF julia> fname = "Test.gpkg" "Test.gpkg"

julia> ds = AG.read(fname) GDAL Dataset (Driver: GPKG/GeoPackage) File(s): Test.gpkg

Number of feature layers: 4 Layer 0: polygons (wkbMultiPolygon) Layer 1: lines (wkbLineString) Layer 2: points (wkbPoint) Layer 3: points_2 (wkbPoint)

julia> dataset_layers = Dict() Dict{Any, Any}()

julia> for i in 0:AG.nlayer(ds)-1 layer = AG.getlayer(ds, i) name_layer = AG.getname(layer) dataset_layers[name_layer] = AG.getlayer(ds, i) end

julia> dataframe = DF.DataFrame(dataset_layers["points"]) 44×3 DataFrame Row │ geom full_id osm_id
│ IGeometr… String String
─────┼───────────────────────────────────────────── 1 │ Geometry: wkbPoint n22062060 22062060 2 │ Geometry: wkbPoint n22062061 22062061 3 │ Geometry: wkbPoint n22062071 22062071 4 │ Geometry: wkbPoint n22062076 22062076 5 │ Geometry: wkbPoint n22062080 22062080 6 │ Geometry: wkbPoint n22066514 22066514 7 │ Geometry: wkbPoint n22066524 22066524 8 │ Geometry: wkbPoint n22066528 22066528 ⋮ │ ⋮ ⋮ ⋮ 38 │ Geometry: wkbPoint n4293533396 4293533396 39 │ Geometry: wkbPoint n4509133611 4509133611 40 │ Geometry: wkbPoint n4965516436 4965516436 41 │ Geometry: wkbPoint n5617933931 5617933931 42 │ Geometry: wkbPoint n6426527642 6426527642 43 │ Geometry: wkbPoint n7921389407 7921389407 44 │ Geometry: wkbPoint n7962854611 7962854611 29 rows omitted

julia> AG.getgeom(dataframe.geom[1], 1) ERROR: GDALError (CE_Failure, code 6): Incompatible geometry for operation

Stacktrace: [1] gdaljl_errorhandler(class::GDAL.CPLErr, errno::Int32, errmsg::Cstring) @ GDAL ~/.julia/packages/GDAL/B6kyy/src/error.jl:36 [2] ogr_g_getgeometryref @ ~/.julia/packages/GDAL/B6kyy/src/GDAL.jl:22147 [inlined] [3] getgeom(geom::ArchGDAL.IGeometry{ArchGDAL.wkbPoint}, i::Int64) @ ArchGDAL ~/.julia/packages/ArchGDAL/zkx2f/src/ogr/geometry.jl:1158 [4] top-level scope @ REPL[8]: Test.zip 1

yeesian commented 3 years ago

Might the following work for you?

using GeoInterface
...
GeoInterface.coordinates(dataframe.geom[1])
pwidas commented 3 years ago

Thanks! Yes this works.