-
Notifications
You must be signed in to change notification settings - Fork 27
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
write geopackage #263
Comments
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 |
Yeah unfortunately so. The implementation in this package should be fixed such that it can write for OGR datasets too |
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. |
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. Works now for GPKG as well :) |
That looks good to me!
How about a name like |
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? |
I think I have a slight preference for updating the write function in Lines 176 to 183 in 04d067e
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 #263 (comment) is still pending and unintuitive. (When |
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. |
The code in your zip file now works under the current master branch. You can use it via |
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
The text was updated successfully, but these errors were encountered: