unexpected behavior for running Proj4.trans on multiple points vs. a single point

[alex-s-gardner]

Hello all... Working with Proj4 and getting some unexpected behavior for running Proj4.trans on multiple points vs. a single point... in my code I've had to put an if, else to switch between single lat lon and multiple lat lon which seems like a bad option. Wondering if anyone has a better way of handling this?

using Proj4

# specify lat and lon values
lat = [69.1, 70.1, 90]
lon = [-49.4, -48, 15]

trans = Proj4.Transformation("EPSG:4326", "EPSG:3413")

# this works [translate multiple lat lon values]
xPt, yPt = trans.(eachrow(hcat(lat, lon)))

# this does not work  [translate single lat lon value]
xPt, yPt = trans.(eachrow(hcat(lat[1], lon[1])))

# this does not work [translate single lat lon value]
xPt, yPt = trans.(hcat(lat[1], lon[1]))

# this works [translate single lat lon value]
xPt, yPt = trans(hcat(lat[1], lon[1]))


# so I have to resort to 
if length(lat) == 1
    xPt, yPt = trans(hcat(lat[1], lon[1]))
  else
    xPt, yPt = trans.(eachrow(hcat(lat, lon)))
  end

[visr]

Hi! It looks like the broadcasting behavior is a bit different than you are expecting here. Without broadcasting, trans takes a point, and returns a point. With broadcasting, it should take multiple points and return multiple points.

So if you look at the output here, it is a list of points:

julia> xPt, yPt = trans.(eachrow(hcat(lat, lon)))
3-element Vector{StaticArrays.SVector{2, Float64}}:
 [-175569.15409751463, -2.281724789444856e6]
 [-113923.50149617658, -2.1737899039931516e6]
 [0.0, 0.0]

Judging from the assigned variables you were expecting an x and a y vector instead, but that is not what is returned.

If your lat and lon coordinate come as separate vectors, and you'd like separate vectors as output as well, you could use such an approach:

# allocate vectors for the output
y = zero(lat)
x = zero(lon)

# fill them one by one
for (i, point) in enumerate(zip(lat, lon))
    x[i], y[i] = trans(point)
end

Now the output looks correct:

julia> x
3-element Vector{Float64}:
 -175569.15409751463
 -113923.50149617658
       0.0

julia> y
3-element Vector{Float64}:
 -2.281724789444856e6
 -2.1737899039931516e6
  0.0

It's probably useful to document this pattern.

[alex-s-gardner]

Ahhh... got it thanks a ton. I'm still dreaming in Matlab so this is really helpful. I also found this discussion on immutable arrays (https://m3g.github.io/JuliaNotes.jl/stable/immutable/) very helpful to grow my understanding of what's going on. In the end I settled on:

specify lat and lon values

lat = [69.1, 70.1, 90]
lon = [-49.4, -48, 15]

trans = Proj4.Transformation("EPSG:4326", "EPSG:3413")
pt = trans.(eachrow(hcat(lat, lon)))

parse x and y values

xPt = first.(pt)
yPt = last.(pt)

[c42f]

xPt = first.(pt)
yPt = last.(pt)

You can also broadcast getindex in these cases (a[i] just means getindex(a, i) in Julia):

xs = getindex.(pt, 1)
ys = getindex.(pt, 2)
zs = getindex.(pt, 3)