write geopackage

[pwidas]

Either I don't understand the documentation or there might be a bug in writing geopackage vector layers to a file. It creates the error: ERROR: LoadError: GDALError (CE_Failure, code 6): Only 1 (Grey/ColorTable), 2 (Grey+Alpha), 3 (RGB) or 4 (RGBA) band dataset supported

Stacktrace: [1] gdaljl_errorhandler(class::GDAL.CPLErr, errno::Int32, errmsg::Cstring) @ GDAL ~/.julia/packages/GDAL/B6kyy/src/error.jl:36 [2] gdalcreatecopy(arg1::Ptr{Nothing}, arg2::String, arg3::Ptr{Nothing}, arg4::Bool, arg5::Ptr{Cstring}, arg6::Base.CFunction, arg7::Ptr{Nothing}) @ GDAL ~/.julia/packages/GDAL/B6kyy/src/GDAL.jl:5902 [3] unsafe_copy(dataset::ArchGDAL.IDataset; filename::String, driver::ArchGDAL.Driver, strict::Bool, options::Ptr{Cstring}, progressfunc::Function, progressdata::Ptr{Nothing}) @ ArchGDAL ~/.julia/packages/ArchGDAL/zkx2f/src/dataset.jl:99 [4] write(dataset::ArchGDAL.IDataset, filename::String; kwargs::Base.Iterators.Pairs{Symbol, ArchGDAL.Driver, Tuple{Symbol}, NamedTuple{(:driver,), Tuple{ArchGDAL.Driver}}}) @ ArchGDAL ~/.julia/packages/ArchGDAL/zkx2f/src/dataset.jl:181 [5] write_ds(fname::String, ds::ArchGDAL.IDataset) @ Main /home/ArchGdal/read_write_gpkg.jl:27 [6] top-level scope @ /home/ArchGdal/read_write_gpkg.jl:34 in expression starting at /home/ArchGdal/read_write_gpkg.jl:34 CPLDestroyMutex: Error = 16 (Device or resource busy)

With the following code, a geopackage dataset is read, and written to the disk (either shapefile or geopackage). Writing to the shapefile works fine. Code and testfile are attached. code_data.zip

[visr]

Here is a slightly reduced example:

julia> import ArchGDAL as AG

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

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


julia> AG.write(ds, "Test_out.gpkg")
ERROR: GDALError (CE_Failure, code 6):
        Only 1 (Grey/ColorTable), 2 (Grey+Alpha), 3 (RGB) or 4 (RGBA) band dataset supported

Stacktrace:
 [1] gdaljl_errorhandler(class::GDAL.CPLErr, errno::Int32, errmsg::Cstring)
   @ GDAL C:\Users\visser_mn\.julia\packages\GDAL\B6kyy\src\error.jl:36
 [2] gdalcreatecopy(arg1::Ptr{Nothing}, arg2::String, arg3::Ptr{Nothing}, arg4::Bool, arg5::Ptr{Cstring}, arg6::Base.CFunction, arg7::Ptr{Nothing})
   @ GDAL C:\Users\visser_mn\.julia\packages\GDAL\B6kyy\src\GDAL.jl:5902
 [3] #unsafe_copy#125
   @ C:\Users\visser_mn\.julia\packages\ArchGDAL\zkx2f\src\dataset.jl:99 [inlined]
 [4] write(dataset::ArchGDAL.IDataset, filename::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ArchGDAL C:\Users\visser_mn\.julia\packages\ArchGDAL\zkx2f\src\dataset.jl:181
 [5] write(dataset::ArchGDAL.IDataset, filename::String)
   @ ArchGDAL C:\Users\visser_mn\.julia\packages\ArchGDAL\zkx2f\src\dataset.jl:181
 [6] top-level scope
   @ REPL[4]:1

The write function for AbstractDataset seems to expect that the dataset is a raster, is that correct?

[yeesian]

The write function for AbstractDataset seems to expect that the dataset is a raster, is that correct?

Yeah unfortunately so. The implementation in this package should be fixed such that it can write for OGR datasets too

[maxfreu]

How could that work? I think the easiest solution would be to detect the dataset driver and then either invoke the old behaviour or direct to a yet-to-write vector writer. The other option would be to have distinct types for raster and vector datasets, but sometimes the boundary between them is blurry, as some tables can contain raster data. So I'd probably go with the first.

The new write functionality could work as follows: create target ds with desired driver, loop over layer definitions in the source ds & recreate each one in the new ds, for each layer loop over all features of the source ds and create the features in the target ds. That would be the wooden hammer.

[maxfreu]

Let me just dump this here quickly:

function AG.deletegeomdefn!(featuredefn::Union{AG.FeatureDefn, AG.IFeatureDefnView}, i::Integer)
    result = GDAL.ogr_fd_deletegeomfielddefn(featuredefn.ptr, i)
    AG.@ogrerr result "Failed to delete geom field $i in the feature definition"
    return featuredefn
end

function copyvectords(ds, name; driver=ArchGDAL.getdriver(ds)) # options
    AG.create(name; driver=driver) do target
        for layeridx in 0:AG.nlayer(ds)-1
            sourcelayer = AG.getlayer(ds, layeridx)
            sourcelayerdef = AG.layerdefn(sourcelayer)
            AG.createlayer(name=AG.getname(sourcelayer),
                           dataset=target,
                           geom=AG.wkbUnknown,  # GPKG only works when using unknown here, deleting the field and recreating it later
                           spatialref=AG.getspatialref(sourcelayer)) do targetlayer
                targetlayerdefn = AG.layerdefn(targetlayer)
                AG.deletegeomdefn!(targetlayerdefn, 0)
                # add geometry field definitions
                for geomfieldidx in 0:AG.ngeom(sourcelayer)-1
                    AG.addgeomdefn!(targetlayer, AG.getgeomdefn(sourcelayerdef, geomfieldidx))
                end
                # add field definitions
                for fieldidx in 0:AG.nfield(sourcelayer)-1
                    AG.addfielddefn!(targetlayer, AG.getfielddefn(sourcelayerdef, fieldidx))
                end
                for feature in sourcelayer
                    AG.addfeature!(targetlayer, feature)
                end
            end
        end
    end
end

AG.create(AG.getdriver("Memory")) do dataset
    layer = AG.createlayer(
        name = "point_out",
        dataset = dataset,
        geom = AG.wkbPoint
    )
    AG.addfielddefn!(layer, "Name", AG.OFTString, nwidth = 32)
    AG.createfeature(layer) do feature
        AG.setfield!(feature, AG.findfieldindex(feature, "Name"), "myname")
        AG.setgeom!(feature, AG.createpoint(100.123, 0.123))
    end
    copyvectords(dataset, "test.sqlite"; driver=AG.getdriver("SQLite"))
end

What do you think? Function naming could be better. ~~Doesn't work yet for gpkg (fails with "GDALError (CE_Failure, code 1): failed to prepare SQL: INSERT INTO "point_out" ( "fid", "", "Name") VALUES (?, ?, ?) - table point_out has no column named"). But works with shp, sqlite.~~

Works now for GPKG as well :)

[yeesian]

That looks good to me!

What do you think? Function naming could be better.

How about a name like writealllayers?

[maxfreu]

Here's the latest version, which works on (at least) shp, gpkg, sqlite and has better performance than the above snippet:

function AG.deletegeomdefn!(featuredefn::Union{AG.FeatureDefn, AG.IFeatureDefnView}, i::Integer)
    result = GDAL.ogr_fd_deletegeomfielddefn(featuredefn.ptr, i)
    AG.@ogrerr result "Failed to delete geom field $i in the feature definition"
    return featuredefn
end

function writevectords(ds, name; driver=ArchGDAL.getdriver(ds), options::Vector{String}=[""], chunksize=20000)
    AG.create(name; driver=driver, options=options) do target
        for layeridx in 0:AG.nlayer(ds)-1
            sourcelayer = AG.getlayer(ds, layeridx)
            sourcelayerdef = AG.layerdefn(sourcelayer)
            ArchGDAL.createlayer(name=AG.getname(sourcelayer),
                                 dataset=target,
                                 geom=AG.wkbUnknown,  # GPKG only works when using unknown here, deleting the field and recreating it later
                                 spatialref=AG.getspatialref(sourcelayer)) do targetlayer
                targetlayerdefn = AG.layerdefn(targetlayer)
                
                # the following try catch block is a bit hacky
                # geopackage files dont work with the default geom field definition,
                # because the driver expects a field "GEOMETRY" but the default is "Geometry"?!
                # on the other hand, the ESRI Shapefile driver does not support GDAL.ogr_l_creategeomfield
                # hence the try catch block as a workaround; it seems to work even after having delted the geomdefn
                try
                    AG.deletegeomdefn!(targetlayerdefn, 0)
                    # add geometry field definitions
                    for geomfieldidx in 0:AG.ngeom(sourcelayer)-1
                        AG.addgeomdefn!(targetlayer, AG.getgeomdefn(sourcelayerdef, geomfieldidx))
                    end
                catch err
                    if !(err isa GDAL.GDALError) && err.msg == "CreateGeomField() not supported by this layer."
                        rethrow(err)
                    end
                end

                # add field definitions
                for fieldidx in 0:AG.nfield(sourcelayer)-1
                    AG.addfielddefn!(targetlayer, AG.getfielddefn(sourcelayerdef, fieldidx))
                end
                
                for chunk in partition(sourcelayer, chunksize)
                    GDAL.ogr_l_starttransaction(targetlayer.ptr)

                    for feature in chunk
                        ArchGDAL.addfeature!(targetlayer, feature)
                    end

                    GDAL.ogr_l_committransaction(targetlayer.ptr)
                end
            end
        end
    end
end

The question is: Where would that code fit best to integrate it?

[yeesian]

The question is: Where would that code fit best to integrate it?

I think writevectordataset(...) in https://github.com/yeesian/ArchGDAL.jl/issues/263#issuecomment-1205438148 might be appropriate in a file like https://github.com/yeesian/ArchGDAL.jl/blob/master/src/dataset.jl

I have a slight preference for updating the write function in https://github.com/yeesian/ArchGDAL.jl/blob/04d067e15b2d9fea2987cb2dd2e130f866dda73c/src/dataset.jl#L176-L183

Something like

target = unsafe_copy(dataset, filename = filename; kwargs...)
for layeridx in 0:AG.nlayer(dataset)-1
    ...
end
destroy(target)

Otherwise the original issue in https://github.com/yeesian/ArchGDAL.jl/issues/263#issue-1045499064 is still pending and unintuitive. (When AG.nlayer(dataset) == 0, it'll still maintain the original behavior.)

[maxfreu]

There are maybe two ways of ipdating the original write function. 1) introduce explicit vector and raster datasets and make a write function for each, or 2) just check the driver with an if else and act according to the driver of the source dataset. I think 1) is overkill and 2) should work fine. The overhead will be negligible.

[maxfreu]

The code in your zip file now works under the current master branch. You can use it via add ArchGDAL#master.