How to create a new shapefile containing a few points

Hi,

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.

2 Likes

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.

1 Like

thank you for your kind suggestions @visr I am checking ArchGDAL.jl, it looks like a very nice and rich library. I was trying to get some inspiration from https://github.com/yeesian/ArchGDAL.jl/blob/master/test/test_cookbook_geometry.jl and GDAL Datasets · ArchGDAL.jl but unfortunately I have not succeed to create my file yet, my skills are limited.

Any help is deeply appreciated. Thanks again!

1 Like

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.

Maybe this (extreme complex ) example can be used as a starting point:

see: ( ‘In [7]:’ ) maybe need adapt …

  • AG.create("pointshapeOut.shp", "ESRI Shapefile") do datasetOut

Can you try the following?

  | | |_| | | | (_| |  |  Version 1.4.2 (2020-05-23)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using ArchGDAL, GDAL

julia> Lon = -44 .+ rand(5); Lat = -23 .+ rand(5)
5-element Array{Float64,1}:
 -22.327553095932217
 -22.179797310271972
 -22.249519172494093
 -22.621065784729858
 -22.364515460011038

julia> ArchGDAL.create(
               "points.shp",
               driver = ArchGDAL.getdriver("ESRI Shapefile")
           ) do ds
           ArchGDAL.createlayer(geom = GDAL.wkbPoint) do layer
               for (lon,lat) in zip(Lon, Lat)
                   ArchGDAL.createfeature(layer) do f
                       ArchGDAL.setgeom!(f, ArchGDAL.createpoint(lat, lon))
                   end
               end
               ArchGDAL.copy(layer, dataset = ds)
           end
       end;

To test that it works:

julia> ArchGDAL.read("points.shp")
GDAL Dataset (Driver: ESRI Shapefile/ESRI Shapefile)
File(s):
  points.shp
  points.shx
  points.dbf

Number of feature layers: 1
  Layer 0: points (wkbPoint)
3 Likes

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!

1 Like

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/>"

Vector drivers — GDAL documentation and ESRI Shapefile / DBF — GDAL documentation has more details. I don’t know if there’s a way to programmatically access that kind of information though :stuck_out_tongue:

4 Likes

Excellent! Thank you again :wave:

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.

Are these codes still working? Do they work if Lon, Lat are Vector{Float64} instead of Array{Float64, 1}? I had trouble reproducing them but I am a very beginner.
I just want to create a shapefile layer from point coordinates.
Thanks anyway

julia> typeof([1.0, 2.0])
Vector{Float64} (alias for Array{Float64, 1})

Thanks for the reply.
Then, I do not understand why the codes do not work when I apply them to 2 Lon, Lat vectors.

1 Like

Probably useful to just create a new thread with a full MWE.

FYI I’m currently working on 2 things that will make this much easier: Shapefile.write and generic geometry wrappers in GeoInterface.jl. Not a huge amount of time available for it though so will be another month or so until its registered.

Thanks, I will

I finally figured out what were the problems. Here is a mininum working example:

# Load the packages
using GDAL, ArchGDAL
# Create Vector{Float64} with 5 point longitudes and latitudes
Lon = [1597315, 1472695, 1499395, 1538365, 1596805]
Lat = [101985.0, 118695.0, -92865.0, 91995.0, 92985.0]

# Check that you can generate spatial points with these coordinates
ArchGDAL.createpoint.(Lon, Lat)

# Use the code posted above to create a "points.shp" shapefile layer:
ArchGDAL.create(
            "points.shp",
            driver = ArchGDAL.getdriver("ESRI Shapefile")
           ) do ds
           ArchGDAL.createlayer(
                    geom = ArchGDAL.wkbPoint,
                    spatialref = ArchGDAL.importEPSG(3978)
                    ) 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 specified that the CRS was EPSG:3978 (appropriate in Québec, Canada).

In my case, the errors in the above code (maybe not errors at the time when posted) came from:

  1. calling geom = GDAL.wkbPoint instead of geom = ArchGDAL.wkbPoint
  2. there was an inversion between lon and lat in createpoint in the initial code posted.

This now works fine.
Thanks

You can also use GeoDataFrames.jl, itself built on top of ArchGDAL. I’m the author.

MWE:

julia> import GeoDataFrames as GDF
julia> import GeoFormatTypes as GFT
julia> using DataFrames
julia> Lon = -44 .+ rand(5); Lat = -23 .+ rand(5)
5-element Vector{Float64}:
 -22.062582685489645
 -22.25804708931404
 -22.13376023130979
 -22.885705294509368
 -22.299512809744762
julia> df = DataFrame(geometry=ArchGDAL.createpoint.(Lon, Lat))
5×1 DataFrame
 Row │ geometry
     │ IGeometr…
─────┼────────────────────
   1 │ Geometry: wkbPoint
   2 │ Geometry: wkbPoint
   3 │ Geometry: wkbPoint
   4 │ Geometry: wkbPoint
   5 │ Geometry: wkbPoint
julia> GDF.write("test.gpkg", df, crs=GFT.EPSG(3978))
"test.gpkg"

edit: Note that we create a DataFrame here, but we essentially accept any table that implements the Tables.jl interface (could be a vector of namedtuples), and the same for any geometry that supports GeoInterface.jl.

Shapefile.jl will write now too.

2 Likes