Failed to create a Shapefile following the GDAL Vector API tutorial

[guillemborrell]

I used an example I found in the tests following the GDAL Vector API tutorial to create a very simple shapefile. The MEMORY driver works fine, but when I try to write the same data down to a file:

ArchGDAL.create("test.shp", driver=ArchGDAL.getdriver("ESRI Shapefile")) do dataset
    layer = ArchGDAL.createlayer(
        name = "point_out",
        dataset = dataset,
        geom = GDAL.wkbPoint
    )
    ArchGDAL.addfielddefn!(layer, "Name", GDAL.OFTString, nwidth = 32)
    featuredefn = ArchGDAL.layerdefn(layer)
    ArchGDAL.createfeature(layer) do feature
        ArchGDAL.setfield!(feature, ArchGDAL.findfieldindex(feature, "Name"), "myname")
        ArchGDAL.setgeom!(feature, ArchGDAL.createpoint(100.123, 0.123))
    end
end

I get the following stack trace (in a jupyter notebook).

Failed to set feature.

Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] macro expansion at /home/guillem/.julia/packages/ArchGDAL/j9NPL/src/utils.jl:14 [inlined]
 [3] setfeature! at /home/guillem/.julia/packages/ArchGDAL/j9NPL/src/ogr/featurelayer.jl:443 [inlined]
 [4] createfeature(::var"#70#72", ::ArchGDAL.IFeatureLayer) at /home/guillem/.julia/packages/ArchGDAL/j9NPL/src/ogr/featurelayer.jl:490
 [5] (::var"#69#71")(::ArchGDAL.Dataset) at ./In[48]:9
 [6] create(::var"#69#71", ::String; kwargs::Base.Iterators.Pairs{Symbol,ArchGDAL.Driver,Tuple{Symbol},NamedTuple{(:driver,),Tuple{ArchGDAL.Driver}}}) at /home/guillem/.julia/packages/ArchGDAL/j9NPL/src/context.jl:216
 [7] top-level scope at In[48]:1
 [8] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091
 [9] execute_code(::String, ::String) at /home/guillem/.julia/packages/IJulia/rWZ9e/src/execute_request.jl:27
 [10] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/guillem/.julia/packages/IJulia/rWZ9e/src/execute_request.jl:86
 [11] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [12] invokelatest at ./essentials.jl:709 [inlined]
 [13] eventloop(::ZMQ.Socket) at /home/guillem/.julia/packages/IJulia/rWZ9e/src/eventloop.jl:8
 [14] (::IJulia.var"#15#18")() at ./task.jl:356

Is this a bug or is it anything I am missing?

Thanks a lot for your work!

[yeesian]

Thanks for giving it a spin! Yeah I can reproduce your example -- unfortunately it looks like a limitation of the shapefile driver in GDAL:

julia> ArchGDAL.create("test.shp", driver=ArchGDAL.getdriver("ESRI Shapefile")) do dataset
          layer = ArchGDAL.createlayer(
              name = "point_out",
              dataset = dataset,
              geom = GDAL.wkbPoint
          )
          ArchGDAL.listcapability(layer)
       end
Dict{String,Bool} with 18 entries:
  "SequentialWrite"    => true
  "DeleteField"        => true
  "IgnoreFields"       => true
  "FastSpatialFilter"  => false
  "DeleteFeature"      => true
  "FastFeatureCount"   => true
  "StringsAsUTF8"      => true
  "CreateGeomField"    => false
  "ReorderFields"      => true
  "MeasuredGeometries" => true
  "FastSetNextByIndex" => true
  "CreateField"        => true
  "RandomWrite"        => true
  "RandomRead"         => true
  "CurveGeometries"    => false
  "FastGetExtent"      => true
  "Transactions"       => false
  "AlterFieldDefn"     => true

One possibility is to create the dataset in memory first, before writing it out as a shapefile, e.g.:

julia> ArchGDAL.create(ArchGDAL.getdriver("MEMORY")) do dataset
           layer = ArchGDAL.createlayer(
               name = "point_out",
               dataset = dataset,
               geom = GDAL.wkbPoint
           )
           ArchGDAL.addfielddefn!(layer, "Name", GDAL.OFTString, nwidth = 32)
           ArchGDAL.createfeature(layer) do feature
               ArchGDAL.setfield!(feature, ArchGDAL.findfieldindex(feature, "Name"), "myname")
               ArchGDAL.setgeom!(feature, ArchGDAL.createpoint(100.123, 0.123))
           end
           ArchGDAL.write(dataset, "test.shp", driver=ArchGDAL.getdriver("ESRI Shapefile"))
       end

[guillemborrell]

Thanks! That's exactly what I needed.

[visr]

This trick of using a MEMORY dataset seems to be needed quite often. I guess in this case it was the "CreateGeomField" capability? Perhaps it's good to reopen this issue as a reminder to document it. On https://yeesian.com/ArchGDAL.jl/dev/features/ there is nothing about feature dataset creation at the moment.

[yeesian]

Yeah agreed, while the design of ArchGDAL is supposed to be permissive with what's possible, it's important to have well-trodden paths (without having to tailor to the particularities of individual drivers) for users who are just looking to get common tasks done.

[ErickChacon]

I just tested the code shared above and I noticed that the geometry is created without name. Is there a way to set up the geometry name? Should the driver do this automatically?

ArchGDAL.create(ArchGDAL.getdriver("MEMORY")) do dataset
   layer = ArchGDAL.createlayer(
       name = "point_out",
       dataset = dataset,
       geom = ArchGDAL.wkbPoint
   )
   ArchGDAL.addfielddefn!(layer, "Name", ArchGDAL.OFTString, nwidth = 32)
   ArchGDAL.createfeature(layer) do feature
       ArchGDAL.setfield!(feature, ArchGDAL.findfieldindex(feature, "Name"), "myname")
       ArchGDAL.setgeom!(feature, ArchGDAL.createpoint(100.123, 0.123))
   end
   ArchGDAL.write(dataset, "test.shp", driver=ArchGDAL.getdriver("ESRI Shapefile"))
end

prueba = AG.read(joinpath("test.shp"))
AG.getlayer(prueba, 0)

# Layer: test
#     Geometry 0 (): [wkbPoint], POINT (100.123 0.123)
#         Field 0 (Name): [OFTString], myname

I also got an error when trying to use GPKG driver. It seems like it is trying to write raster data. Any suggestion?

ArchGDAL.create(ArchGDAL.getdriver("MEMORY")) do dataset
   layer = ArchGDAL.createlayer(
       name = "point_out",
       dataset = dataset,
       geom = ArchGDAL.wkbPoint
   )
   ArchGDAL.addfielddefn!(layer, "Name", ArchGDAL.OFTString, nwidth = 32)
   ArchGDAL.createfeature(layer) do feature
       ArchGDAL.setfield!(feature, ArchGDAL.findfieldindex(feature, "Name"), "myname")
       ArchGDAL.setgeom!(feature, ArchGDAL.createpoint(100.123, 0.123))
   end
   ArchGDAL.write(dataset, "test.gpkg", driver=ArchGDAL.getdriver("GPKG"))
end

# ERROR: GDALError (CE_Failure, code 6):
#     Only 1 (Grey/ColorTable), 2 (Grey+Alpha), 3 (RGB) or 4 (RGBA) band dataset supported