Extracting coordinates from Shape File

Specific Domains Geo
,
1

Hi Community

I tried to read in a Shape File with Austrian Boarders typing

julia> A=Shapefile.Table(Path_Shapefile)
Shapefile.Table{Union{Missing, Shapefile.Polygon}} with 1 rows and the following 16 columns:

geometry, ID, LANGUAGE, NAME, PARENT_ID, ADMN_LEVEL, ADM_LEV_DE, ADM_LEV_EN, ADM_LEV, MAPS_LEVEL, MapLev_min, MapLev_max, ID_ORIG, FEATTYP, TZID, UTC_OFFSET

julia> B=A.geometry
1-element Vector{Union{Missing, Shapefile.Polygon}}:
 Polygon(5890 Points)

How is it now possible to get the Polygon values stored in variable B ?

2

You accessed the geometry column of the table. You can access any other column with similar syntax. I am assuming that is how Shapefile.Table works. If it doesn´t work that way, you can also try GeoTables.jl to load the shapefile with pure Julia geometries:

using GeoTables

table = GeoTables.load("myfile.shp")

table.geometry
table.ADMN_LEVEL
...

Notice that GeoTables.jl will recognize any file format, not just Shapefile. It works with Geopackage or any other file format supported by GDAL.

1 Like
3

The data types in Shapefile.jl are defined here: https://github.com/JuliaGeo/Shapefile.jl/blob/master/src/Shapefile.jl#L58 , so B.points gives you a vector of Points while B.parts denotes the indices at which a new polygon starts.

4

Thanks fabiangans and juliohm,

GeoTables is able to open the shapefile

julia> using GeoTables

julia> table = GeoTables.load(Path_Shapefile)
1 GeoTable{2,Float64}
  variables
    └─ADMN_LEVEL (Int64)
    └─ADM_LEV (String)
    └─ADM_LEV_DE (String)
    └─ADM_LEV_EN (String)
    └─FEATTYP (Int64)
    └─ID (String)
    └─ID_ORIG (String)
    └─LANGUAGE (String)
    └─MAPS_LEVEL (Int64)
    └─MapLev_max (Int64)
    └─MapLev_min (Int64)
    └─NAME (String)
    └─PARENT_ID (String)
    └─TZID (String)
    └─UTC_OFFSET (String)
  domain: 1 GeometrySet{2,Float64}

julia> Geo = table.table.geometry
1-element Vector{Union{Missing, Shapefile.Polygon}}:
 Polygon(5890 Points)

But here is the same question how to receive the coordinates of the 5980 points stored in Geo.

5

Hi fabiangans,

Unfortunately I receive this error message when trying your suggestion:

julia> B
1-element Vector{Union{Missing, Shapefile.Polygon}}:
 Polygon(5890 Points)

julia> B.points
ERROR: type Array has no field points
Stacktrace:
 [1] getproperty(x::Vector{Union{Missing, Shapefile.Polygon}}, f::Symbol)
   @ Base ./Base.jl:33
 [2] top-level scope
   @ none:1
6

@Dieter take a look at the Meshes.jl documentation which is where all these geometries are defined. You can collect the vector of geometries, pick any geometry in the vector and collect its vertices:

# column of geometries
geoms = table.geometry

# first geometry
g = geoms[1]

# list of vertices
vs = vertices(g)

# list of coordinates
xs = coordinates.(vs)
7

Sorry for the error, B.points is a vector of length one, so B[1].points should access the points field of the first element.

8

Thanks @fabiangans, now it works!

9

Update for those coming late here:

Please replace GeoTables.load by GeoIO.load.

GeoIO.jl is the new home of the load and save functions.