I have the latitude and longitude of a few points (in WGS84, EPS 4326) and I would like to save my points in a Esri shapefile.
My 5 points have coordinates like Lon = -44 .+ rand(5); Lat = -23 .+ rand(5)
Please, how could I do this in Julia?
Thank you very much for your attention.
PS. I was checking Shapefile.jl and found I can read a shape files using
table = Shapefile.Table(myFileName.shp)
but unfortunately I could not find and simple example of how to create a file with my data.
Shapefile.jl is so far for reading only, writing is not supported. If GeoJSON is also an acceptable file format, you can write with GeoJSON.jl. For writing Shapefiles you can use ArchGDAL.jl. There it is supported, though I don’t have a simple example ready.
Could you share what you have so far in ArchGDAL.jl? Maybe we can help make it work. These kind of examples would be great to have in the documentation.
The above code is almost perfect: it creates 3 out of 4 files I would expect: points.shp points.dbf and points.shx (only the projection .prj file is missing). Please, do you know where I could try to set spatialref = ArchGDAL.importEPSG(4326) to see if I get the missing .prj file?
As a minor note, I replaced ArchGDAL.createpoint(lat, lon) with ArchGDAL.createpoint(lon,lat), thus the points are the right geographical location (setting manually the projection code 4326 using an external software).
Thank you again guys for all your answers and effort!
I tried as well but could not manage. Putting this at the end before closing the dataset gave the following error:
spatialref = ArchGDAL.importEPSG(4326)
wkt_crs = ArchGDAL.toWKT(spatialref)
ArchGDAL.setproj!(ds, wkt_crs)
# GDAL error: points: Dataset does not support the SetSpatialRef() method
Ah, you’ll have to specify the spatial reference when defining the layer:
ArchGDAL.create(
"points.shp",
driver = ArchGDAL.getdriver("ESRI Shapefile")
) do ds
ArchGDAL.createlayer(
geom = GDAL.wkbPoint,
spatialref = ArchGDAL.importEPSG(4326)
) do layer
for (lon,lat) in zip(Lon, Lat)
ArchGDAL.createfeature(layer) do f
ArchGDAL.setgeom!(f, ArchGDAL.createpoint(lon, lat))
end
end
ArchGDAL.copy(layer, dataset = ds)
end
end;
I think some of these capabilities are driver-specific, and I had to experiment to find out too. See e.g.
julia> ArchGDAL.create(
ArchGDAL.getdriver("ESRI Shapefile")
) do ds
ArchGDAL.listcapability(ds)
end
Dict{String,Bool} with 6 entries:
"CreateLayer" => true
"DeleteLayer" => true
"CreateGeomFieldAfterCreateLayer" => false
"CurveGeometries" => false
"Transactions" => false
"EmulatedTransactions" => false
julia> ArchGDAL.options(
options(drv::ArchGDAL.Driver) in ArchGDAL at /Users/yeesian/.julia/packages/ArchGDAL/ohMpc/src/driver.jl:48
julia> ArchGDAL.options(ArchGDAL.getdriver("ESRI Shapefile"))
"<CreationOptionList/>"
Apologies for resurrecting an old thread, but just wanted to add that these instructions worked for shapefile output, but gave an error with GeoPackage output; I had to replace ArchGDAL.createfeature() with ArchGDAL.addfeature() to get it to work for GeoPackage.