evetion / wellknowngeometry.jl Goto Github PK

View Code? Open in Web Editor NEW

3.0 4.0 2.0 299 KB

Reading and writing Well Known Text (WKT) and Well Known Binary (WKB)

Home Page: https://www.evetion.nl/WellKnownGeometry.jl/dev

License: MIT License

Julia 100.00%

gis julia wellknowntext wellknownbinary wellknowngeometry io wkb wkt

wellknowngeometry.jl's Introduction

WellKnownGeometry

Reading and writing Well Known Text (WKT) and Well Known Binary (WKB) based on GeoInterface.jl. See this Wikipedia page for an explanation of Well-known text and binary geometry.

Given a GeoInterface compatible geometry, this package can generate the WKT and WKB representation of it. It also does the reverse, as it implements GeoInterface for WKT or WKB strings.

WellKnownText

Given a WKT string, we can retrieve the type and underlying coordinates, and thus convert it to other geometries using GeoInterface. Note that WKT strings are wrapped by GeoFormatTypes so we can distinguish them from any other strings.

using ArchGDAL
using GeoFormatTypes
using WellKnownGeometry
using GeoInterface

wkts = "POINT (30 10)"
wkt = GeoFormatTypes.WellKnownText(GeoFormatTypes.Geom(), wkts)

GeoInterface.geomtrait(wkt)  # PointTrait()
GeoInterface.ncoord(wkt)  # 2
GeoInterface.coordinates(wkt)  # 2-element Vector{Float64}: 30.0 10.0

p = convert(ArchGDAL.IGeometry{ArchGDAL.wkbPoint}, wkt)  # Geometry: POINT (30 10)

As ArchGDAL geometries implement GeoInterface, we can generate the WKT for it.

wkt = WellKnownGeometry.getwkt(p)  # WellKnownText{GeoFormatTypes.Geom}(GeoFormatTypes.Geom(), "POINT (30.0 10.0)")
GeoFormatTypes.val(wkt)  # "POINT (30.0 10.0)"

getwkt

WellKnownBinary

Given a WKB byte string, we can retrieve the type and underlying coordinates, and thus convert it to other geometries using GeoInterface. Note that WKB byte strings are wrapped by GeoFormatTypes so we can distinguish them from any other byte strings.

using ArchGDAL
using GeoFormatTypes
using WellKnownGeometry
using GeoInterface

wkbs = [0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40]
wkb = GeoFormatTypes.WellKnownBinary(GeoFormatTypes.Geom(), wkbs)

GeoInterface.geomtrait(wkb)  # PointTrait()
GeoInterface.ncoord(wkb)  # 2
GeoInterface.coordinates(wkb)  # 2-element Vector{Float64}: 30.0 10.0

p = convert(ArchGDAL.IGeometry{ArchGDAL.wkbPoint}, wkb)  # Geometry: POINT (30 10)

As ArchGDAL geometries implement GeoInterface, we can generate the WKB for it.

wkb = WellKnownGeometry.getwkb(p)  # WellKnownBinary{GeoFormatTypes.Geom, Vector{UInt8}}(GeoFormatTypes.Geom(), UInt8[0x01, ..., 0x40])
GeoFormatTypes.val(wkb)  # 21-element Vector{UInt8}:  0x01 0x01 ... 0x00 0x40

getwkb

wellknowngeometry.jl's People

Contributors

Stargazers

Watchers

Forkers

visr asinghvi17

wellknowngeometry.jl's Issues

Support empty geometries

Currently, we support EMPTY in WKT, but ArchGDAL returns [0.0, 0.0], while we return []. For WKB, it's worse, as there is no fixed definition on how to represent WKB. See https://blog.cleverelephant.ca/2010/03/nothing-nada-zip-bupkus.html

Parser does not work with views

MWE:

julia> import GeoFormatTypes as GFT, WellKnownGeometry as WKG, GeoInterface as GI

julia> geoms = [UInt8[0x47, 0x50, 0x00, 0x01, 0x84, 0x08, 0x00, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0xb0, 0x41, 0xd2, 0x8e, 0xc5, 0xa9, 0x15, 0x41, 0x40, 0xdb, 0xcf, 0xb9, 0xb0, 0x82, 0x08, 0x41]]

julia> GFT.WellKnownBinary(GFT.Geom(), view(geoms[1], (8+1):length(geoms[1])))
GeoFormatTypes.WellKnownBinary{GeoFormatTypes.Geom, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{Int64}}, true}}(GeoFormatTypes.Geom(), UInt8[0x01, 0x01, 0x00, 0x00, 0x00, 0xb0, 0x41, 0xd2, 0x8e, 0xc5  …  0x15, 0x41, 0x40, 0xdb, 0xcf, 0xb9, 0xb0, 0x82, 0x08, 0x41])

julia> GI.trait(ans)

julia> GFT.WellKnownBinary(GFT.Geom(), getindex(geoms[1], (8+1):length(geoms[1])))
GeoFormatTypes.WellKnownBinary{GeoFormatTypes.Geom, Vector{UInt8}}(GeoFormatTypes.Geom(), UInt8[0x01, 0x01, 0x00, 0x00, 0x00, 0xb0, 0x41, 0xd2, 0x8e, 0xc5  …  0x15, 0x41, 0x40, 0xdb, 0xcf, 0xb9, 0xb0, 0x82, 0x08, 0x41])

julia> GI.trait(ans)
PointTrait()

Non-Float64 types don't work correctly

As their byte representation is shorter.

Backport Findall to support Julia 1.6

Any reason to require julia 1.7 instead of 1.6? a7c0250
This could become a pretty important package to the ecosystem, which mostly still supports the LTS release.
And what do you think about moving it to JuliaGeo?

Originally posted by @visr in #9 (comment)

Support Extended variants

Those with a CRS prefixed, both in Binary and Text form.

Parametrized parser for speed

It would be great to have a parameterized parser object for speed. Things we could specialize on are:

Endianness (little->big endian is just bswap on the result Float64 or Int)
Geometry type: when parsing SQL query results, we often know that the column contains only one kind of geometry. In this case, it would be nice to be able to call a parser which only parses, e.g., points, or polygons, etc. This parameterized parser could simply error if the geometry was found to be of an incorrect type. This also guarantees return type and would probably make fast inner loops possible when parsing columns in GeoParquet and GeoPackage.
- The generic parser would be type unstable and do all necessary type checking, similar to the current approach.

Efficient coordinate collection from `LineString`s using `reinterpret`

import GeoFormatTypes as GFT, WellKnownGeometry as WKG, GeoInterface as GI

tups = tuple.(rand(300_000), rand(300_000))
geoms_as_points = GFT.val.(WKG.getwkb.(GI.Point.(tups)))
geoms_as_linestring = GFT.val(WKG.getwkb(GI.LineString(tups)))

# parse geoms as points, each WKB individually
@benchmark begin
    map($(geoms_as_points)) do geom
        only(reinterpret(Tuple{Float64, Float64}, view(geom, (5+1):length(geom)))) # this has to be parsed differently in EWKB and GeoPkg wkb
    end
end # median time 3ms v/s regular parsing 300ms

# parse the whole linestring
@benchmark collect(reinterpret(Tuple{Float64, Float64}, view($ls_as_wkb, (1+4+4+1):length($ls_as_wkb))))
# median time 290ns vs unknown (still running) 
# we could even drop collect here

Support ZM geometries

Now we only support 2d (XY) ones.

Broken doc links

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

`getwkb` should return a WellKnownGeometry

Issue parsing geometry when no space exists between type and data

Trying to parse this WKT file line-by-line using WellKnownGeometry/GeoFormatTypes fails, returning a missing geometry.
small_example.wkt.txt

In the file, the syntax goes MULTIPOLYGON(...) with no space between MULTIPOLYGON and (...). This seems to throw the parser off.

A minimal test to replicate this is the following:

using WellKnownGeometry, GeoFormatTypes, GeoInterface

julia> wkts = "POINT (30 10)"
"POINT (30 10)"

julia> wkt = GeoFormatTypes.WellKnownText(GeoFormatTypes.Geom(), wkts)
WellKnownText{GeoFormatTypes.Geom}(GeoFormatTypes.Geom(), "POINT (30 10)")

julia> GeoInterface.geomtrait(wkt)  # PointTrait()
PointTrait()

julia> GeoInterface.ncoord(wkt)  # 2
2

julia> GeoInterface.coordinates(wkt)  # 2-element Vector{Float64}: 30.0 10.0
2-element Vector{Float64}:
 30.0
 10.0

julia> wkts = "POINT(30 10)"
"POINT(30 10)"

julia> wkt = GeoFormatTypes.WellKnownText(GeoFormatTypes.Geom(), wkts)
WellKnownText{GeoFormatTypes.Geom}(GeoFormatTypes.Geom(), "POINT(30 10)")

julia> GeoInterface.geomtrait(wkt)
┌ Warning: unknown geometry type
│   geom.val = "POINT(30 10)"
└ @ WellKnownGeometry ~/.julia/packages/WellKnownGeometry/t1n94/src/wkt.jl:112

Ideal behaviour would be to understand and parse this correctly. LibGEOS seems to do this in the way I would expect (and as outlined in this issue) so I am using that for now, but would like to switch to a Julia solution.

GeoInterface.getcoord error on LineString

I stumbled on this error:

import GeoFormatTypes as GFT
using GeoInterface
wkt = GFT.WellKnownText(GFT.Geom(), "POINT (30 10)")
GeoInterface.coordinates(wkt)  # works
wkt = GFT.WellKnownText(GFT.Geom(), "LINESTRING (30.0 10.0, 10.0 30.0, 40.0 40.0)")
GeoInterface.coordinates(wkt)  # fails

ERROR: ArgumentError: cannot parse " " as Float64
Stacktrace:
  [1] _parse_failure(T::Type, s::String, startpos::Int64, endpos::Int64) (repeats 2 times)
    @ Base .\parse.jl:373
  [2] #tryparse_internal#477
    @ .\parse.jl:369 [inlined]
  [3] tryparse_internal
    @ .\parse.jl:366 [inlined]
  [4] #parse#478
    @ .\parse.jl:379 [inlined]
  [5] parse
    @ .\parse.jl:379 [inlined]
  [6] getcoord(#unused#::PointTrait, geom::GeoFormatTypes.WellKnownText{GeoFormatTypes.Geom}, i::Int64)
    @ WellKnownGeometry d:\visser_mn\.julia\dev\WellKnownGeometry\src\wkt.jl:132

Stripping it down a bit, it seems to stumble on the comma or space before the second point in the linestring

wkt = GFT.WellKnownText(GFT.Geom(), "LINESTRING (30.0 10.0, 10.0 30.0, 40.0 40.0)")
p = GeoInterface.getgeom(wkt, 2)
GeoInterface.getcoord(p, 1)  # parse error
GeoInterface.getcoord(p, 2)  # 10 (should be 30)

Add `GeoInterface.convert` methods

Would be nice if convert worked ;)

We can add a convert(WellKnownGeometry, geom) method now too.