acgold / SimpleFeatures.jl

Working with simple feature GIS data in Julia
https://acgold.github.io/SimpleFeatures.jl/
MIT License
9 stars 0 forks source link

st_read not working for SpatiaLite files #29

Closed GerritGue closed 1 year ago

GerritGue commented 1 year ago

Hi,

first of all, thank you for your efforts to bring the functionality of the sf package to Julia!

Describe the bug

While diving into the package, I ran into an error regarding SpatiaLite (".sqlite") files. I just created a dummy file containing random points which can be downloaded here.

Minimal Working Example

import SimpleFeatures as sf

random_points = sf.st_read("random_points.sqlite")

I get the following error message:

ERROR: ArgumentError: column name :geom not found in the data frame
Stacktrace:
  [1] lookupname
    @ ~/.julia/packages/DataFrames/dgZn3/src/other/index.jl:396 [inlined]
  [2] getindex
    @ ~/.julia/packages/DataFrames/dgZn3/src/other/index.jl:405 [inlined]
  [3] parentcols
    @ ~/.julia/packages/DataFrames/dgZn3/src/other/index.jl:455 [inlined]
  [4] getindex
    @ ~/.julia/packages/DataFrames/dgZn3/src/dataframerow/dataframerow.jl:212 [inlined]
  [5] getproperty(r::DataFrameRow{DataFrame, DataFrames.Index}, idx::Symbol)
    @ DataFrames ~/.julia/packages/DataFrames/dgZn3/src/dataframerow/dataframerow.jl:302
  [6] (::SimpleFeatures.var"#19#20"{ArchGDAL.Dataset})(table::ArchGDAL.FeatureLayer)
    @ SimpleFeatures ~/.julia/packages/SimpleFeatures/eNsu2/src/from_gdf.jl:105
  [7] getlayer(::SimpleFeatures.var"#19#20"{ArchGDAL.Dataset}, ::ArchGDAL.Dataset, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArchGDAL ~/.julia/packages/ArchGDAL/AZONk/src/context.jl:267
  [8] getlayer
    @ ~/.julia/packages/ArchGDAL/AZONk/src/context.jl:264 [inlined]
  [9] st_read
    @ ~/.julia/packages/SimpleFeatures/eNsu2/src/from_gdf.jl:91 [inlined]
 [10] (::SimpleFeatures.var"#14#15")(ds::ArchGDAL.Dataset)
    @ SimpleFeatures ~/.julia/packages/SimpleFeatures/eNsu2/src/from_gdf.jl:78
 [11] read(f::SimpleFeatures.var"#14#15", args::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArchGDAL ~/.julia/packages/ArchGDAL/AZONk/src/context.jl:267
 [12] read
    @ ~/.julia/packages/ArchGDAL/AZONk/src/context.jl:264 [inlined]
 [13] #st_read#13
    @ ~/.julia/packages/SimpleFeatures/eNsu2/src/from_gdf.jl:73 [inlined]
 [14] st_read(fn::String)
    @ SimpleFeatures ~/.julia/packages/SimpleFeatures/eNsu2/src/from_gdf.jl:72

Reading the data with ArchGDAL works as expected

using ArchGDAL
using DataFrames

random_points = ArchGDAL.read("random_points.sqlite")

DataFrames.DataFrame(ArchGDAL.getlayer(random_points, 0))

leads to

100×2 DataFrame
 Row │                     id    
     │ IGeometr…           Int32 
─────┼───────────────────────────
   1 │ Geometry: wkbPoint      0
   2 │ Geometry: wkbPoint      1
   3 │ Geometry: wkbPoint      2
   4 │ Geometry: wkbPoint      3
   5 │ Geometry: wkbPoint      4
   6 │ Geometry: wkbPoint      5
   7 │ Geometry: wkbPoint      6

Reading GeoJSON and GeoPackage files works as expected.

Edit: I also tried to load the file using GeoDataFrames and it works as expected:

import GeoDataFrames as GDF

GDF.read("random_points.sqlite")

leads to

100×2 DataFrame
 Row │ GEOMETRY            id    
     │ IGeometr…           Int64 
─────┼───────────────────────────
   1 │ Geometry: wkbPoint      0
   2 │ Geometry: wkbPoint      1
   3 │ Geometry: wkbPoint      2
   4 │ Geometry: wkbPoint      3
   5 │ Geometry: wkbPoint      4
   6 │ Geometry: wkbPoint      5

Versioninfo()

Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 7 3700X 8-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
  Threads: 10 on 16 virtual cores
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 10

And I am using SimpleFeatures.jl v0.2.2

acgold commented 1 year ago

Hi,

Thanks for opening an issue. The geometry column in SimpleFeatures should generally be called "geom". Issues arise if this is not the case because in some instances I (lazily) declared the geometry column to always be named "geom" rather than providing flexibility. Maybe that is the issue? One workaround is to read the file with ArchGDAL and convert it to a SimpleFeature object (https://acgold.github.io/SimpleFeatures.jl/dev/reference/#SimpleFeatures.df_to_sf).

I created SimpleFeatures so I could test out Julia for a project focused on a large commuting dataset, so there are some cases where I prioritized functionality for my use case rather than general use. I have since left my postdoc position and gotten a new job, and I am not able to keep this package up to date as the Julia GIS ecosystem is maturing. PRs are very welcome, but I would also highly recommend taking a look at GeoDataframes and ArchGDAL which have much better long-term support. It looks like GeoDataframes recently added the ability to hold CRS data within DataFrame metadata (https://github.com/evetion/GeoDataFrames.jl/releases/tag/v0.3.2), which was the functionality that was previously lacking and spurred me to create SimpleFeatures.

GerritGue commented 1 year ago

Thank your for your reply, this solves the issue and congrats on your new position