CoordinateTransformations.jl API

[visr]

Not yet finished, but parking the work in progress PR here until I have time to finish it. Part of #44, see https://github.com/JuliaGeo/Proj4.jl/pull/44#issuecomment-759693533

TODO

• [x] documentation (input can be any indexable length 2-4, output is same length, but SVector{Float64}), network demo (quickstart in README done)
• [x] Base.show(::Transformation)
• [x] Add methods for vectors of points, use proj_trans_array / proj_trans_generic? (better in a separate PR, if at all)
• [x] Implement CoordinateTransformations.compose

[visr]

I'm looking for some input on input, related to this bit of code: https://github.com/JuliaGeo/Proj4.jl/blob/366a5ff4046bcf73564ae51815d51160d868bd30/src/coord.jl#L74-L112

The PROJ C API always works with length 4 arrays, where the first two positions are x or y, then height and time. For when you only care about the x and y, we default to a height of 0 and infinite time, just like the cs2cs application. In cs2cs, if your input is size 4, you get size 4 back. If the input is size 2 or 3, you get size 3 back, always including height. My idea was to instead always pass back the same size as your input. If you care about height, you will only get it back if you pass it in as well. Does that sound reasonable? Or is it better to follow cs2cs here, and always return height, such that users are more aware if the operation affects the height as well.

Besides the length, there is also the matter what types to accept and return. Right now, I always return SVector{Float64}, with length as above. If the input is an SVector as well, I can specialize on the length, otherwise, for e.g. AbstractVector/Tuple input I check the length and adjust the output length accordingly. Though I guess that would mean it is type unstable. Is this the right approach? Any other ideas?

Alternatively I could try to give back the same type that comes in, but with eltype Float64. E.g. people passing tuples might prefer a tuple back, and people passing a GeometryBasics.Point might prefer that back, using a method like f(x) = typeof(x)(SA_F64[x...]) which can return that. Or would this make things too complicated?

Would be happy to get ideas from anyone, tagging @c42f since he helped out before.

[c42f]

My idea was to instead always pass back the same size as your input. If you care about height, you will only get it back if you pass it in as well. Does that sound reasonable?

Sounds reasonable to me. As an extension of this, you could check that any dummy time and height values passed are returned back unchanged. If not, throw an error. (Optionally have a flag in the transformation to turn off such checking if performance concerns warrant it?)

If the input is an SVector as well, I can specialize on the length, otherwise, for e.g. AbstractVector/Tuple input I check the length and adjust the output length accordingly. Though I guess that would mean it is type unstable. Is this the right approach? Any other ideas?

I think the right approach here is to use StaticArrays.similar_type for any input which is a subtype of StaticVector (ie, generalize your functions which currently take only SVector)

For other AbstractVector, either:

• Type instability may not be so bad if the compiler can figure out that the return type is Union{SVector{2,Float64},SVector{3,Float64},SVector{4,Float64}}? That way the small union optimization may kick in. The compiler can probably figure this out but you could add a type assertion if necessary. Do a performance test?
• Failing that, you should probably just use Base.similar() to construct an appropriate array which can be filled.

[visr]

Thanks for the helpful tips! I tried to implement the changes in d60567c17d722b4fe8fcbd504286c9deaac7604a.

As an extension of this, you could check that any dummy time and height values passed are returned back unchanged. If not, throw an error.

I'm afraid this would get in the way for 2D applications that don't mind also if the height is changed. I believe the time always comes back unchanged.

Type instability may not be so bad if the compiler can figure out that the return type is Union{SVector{2,Float64},SVector{3,Float64},SVector{4,Float64}}?

For a size 2 Vector{Float64} input, running with Test.@inferred I get

SVector{2, Float64} does not match inferred return type SArray{S, Float64, 1, L} where {S<:Tuple, L}. Since SVector{2, Float64} is an alias for SArray{Tuple{2}, Float64, 1, 2}), I guess this means the compiler did not figure out that it is 2, 3, or 4, but any n? Not sure if that is worse than a Union of 3 types.

BenchmarkTools shows it is almost equally fast (1.3 μs) for SVector, Vector and Tuple input. Only for Vector it gives a single allocation. So I guess that is good enough.

[c42f]

I'm afraid this would get in the way for 2D applications that don't mind also if the height is changed.

True. In that case I guess what you've got is fine.

I believe the time always comes back unchanged.

Ah of course... it would be surprising if proj supported nontrivial temporal reference frames / datum shifts :laughing:

SVector{2, Float64} does not match inferred return type SArray{S, Float64, 1, L} where {S<:Tuple, L}. Since SVector{2, Float64} is an alias for SArray{Tuple{2}, Float64, 1, 2}), I guess this means the compiler did not figure out that it is 2, 3, or 4, but any n? Not sure if that is worse than a Union of 3 types.

Ok, good to know. In that case, I'd suggest adding a type assertion with the Union of three types and checking whether that helps. You'll also want to benchmark actually doing something with the returned value (eg, something simple such as accumulating it with +), because that's where the performance hit will be. I think a union of three types should be better than the UnionAll which the compiler is inferring.

[visr]

I'd suggest adding a type assertion with the Union of three types and checking whether that helps

It doesn't seem to help. I updated the benchmark, and the performance of vectors is really quite close, so perhaps it's fine to leave it. This is the benchmark I used:

using StaticArrays
using BenchmarkTools
using Proj4
using Test

const trans = Proj4.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)

# test different input types
ps = SA[5.39, 52.16]
pv = [5.39, 52.16]
pt = (5.39, 52.16)

function timer(p)
    tot = 0.0
    for _ in 1:10_000
        tot += sum(trans(p))
    end
    tot
end

@btime timer($ps) # 13.512 ms (0 allocations: 0 bytes)
@btime timer($pv) # 16.543 ms (30000 allocations: 625.00 KiB)
@btime timer($pt) # 13.179 ms (0 allocations: 0 bytes)

Test.@inferred trans(ps)  # pass
Test.@inferred trans(pv)  # return type SVector{2, Float64} does not match inferred return type SArray{S, Float64, 1, L} where {S<:Tuple, L}
Test.@inferred trans(pt)  # pass

Overall I think this PR is ready. Does anyone have other feedback or is this good to merge? Will merge and tag in a few days without comments.

[c42f]

It doesn't seem to help. I updated the benchmark, and the performance of vectors is really quite close, so perhaps it's fine to leave it. This is the benchmark I used:

In that case looks good to me. It seems that the cost of the transformation is just relatively large compared to other costs in this code.

I also tried out Geodesy.UTMfromLLA and it was slightly slower than the comparable transformation from Proj so it seems that the cost of interop with the proj C library is well under control here. (Albeit being slightly disappointing in other aspects :sweat_smile: )

[c42f]

Nice!