Hello, I am fairly new to Julia and absolutely new to ArchGDAL, so please excuse me if what I am reporting is actually a misconception on my part. I would like to first thank the developers of both the language and the package for their efforts.
I am not, however, new to programming or spatial analysis and I have implemented in R, python and C++ uses of GDAL, and never found this situation. I would like to know if this is a bug or there is something I do not understand.
The issue is that as a general rule in the other languages/packages, the x coordinate even when using EPSG:4326 is the lon coordinate and y is the lat. This allows to generalize functions, and in ArchGDAL it seems to not be the case. Let me explain myself better and provide an example.
I use the following libraries:
using PyPlot
using ArchGDAL
const AG = ArchGDAL
For convenience, and because I am not a fan of the do…end block syntax, I created a function to reproject a set of points, that I might be taking out of a CSV or whatever other source:
function reproject_points(x,y,EPSG_in,EPSG_out)
if length(size(x))>0
source = AG.importEPSG(EPSG_in)
target = AG.importEPSG(EPSG_out)
point = AG.createpoint.(x,y)
AG.createcoordtrans(source, target) do transform
AG.transform!.(point, Ref(transform))
end
else
source = AG.importEPSG(EPSG_in)
target = AG.importEPSG(EPSG_out)
point = AG.createpoint(x,y)
AG.createcoordtrans(source, target) do transform
AG.transform!(point, transform)
end
end
return point
end
This is the data I will be using for this example, it is point coordinates in EPSG:2062
xs = [607206.9482826635,600884.2090243029,601285.019884102,601271.2564008039,599963.1503200289,599664.2732100531,601222.7195474328,605584.0543730418,601868.9331501597,600657.6185030371,599676.4642521591,599710.787271369,601785.6991023883,600381.0672790125,605614.9847654601,605676.685498099,604895.3714377616,607844.9892708402,597517.9261878692,600189.9270375178,608901.2223834917,606680.0393119903,605875.0781300119,600990.2953122013,599709.8854868729,600071.8545435866,597960.0596096516,599140.8002564722,598947.6141444272,598619.3702475202,599036.2309225708,598305.7552071137,599852.1820001382,599385.9559182358,599420.008648855,606696.4380820625,606713.4575155701,606750.1489200754,606770.9789345554,606912.3133015297,597999.1989373707,597444.1489742319,599243.5713090095,600985.4180041264,597821.5931688325,597614.1835169938,598492.4293076693,600192.4704526738,600881.9574123736,598981.8396000165,599061.6798259569,599064.2043957161,600664.0548508466,599099.4885307772,598282.9211607504,600015.818270395,599555.781922397,599710.6595791673,599863.521450467,599350.4881577014,599755.9993486871,597602.4557140957,599108.1827850383,599102.5838569351,599084.966560575,605271.9340692578,598664.3126318238,598956.7220973949,598611.2321848627,598339.7573898082,598934.5454809486,598831.9463030831,599838.0321587456,598625.2015675211,597595.6818008538,600997.920543785,598919.3913176566,597767.64002293,599303.402291598,599095.2671827096,599594.5563932064,597702.4362196913,598282.4404220523,598971.595085984,598401.9905808555,599236.127706109,599306.6780913315,599082.5470946438,599434.6384900511,598243.2548718995,599424.9037081582,599578.0048930706,599559.4238164076,599691.9632582335,600071.9750061606,598994.7992599889,599262.0813262118,599075.9069290818,598643.908856648,597181.7616042689,598013.4985272097]
ys = [646790.5658288356,648509.1337229714,648367.3547832256,648338.270249308,647448.7987419742,646562.1648201225,648668.6353989807,642167.5543795546,647681.4240737665,646936.5952231894,637786.1488118544,639416.1262479564,647590.8466447088,646776.56734984,641530.2915909081,641691.6452882218,647596.0027563099,649298.5971032989,638065.1676245964,646684.1615885856,652515.3428531999,640778.1341588998,642265.7398512278,638068.1269821333,639016.4351534036,639372.9674865914,638391.3461249629,639727.4583385806,638739.7379113979,638545.3399807084,639294.9873405725,638525.9171089388,639609.0731436828,639941.195687185,639672.6726679131,642656.7993175152,642524.0348008994,642560.2917718817,642451.5382035525,642137.8804459467,638369.6676396746,638043.6111432577,639102.7515251131,638438.7964159701,638207.4155945882,638408.643140416,638807.1455215111,638121.893765639,638507.937714509,639350.5204853417,639492.0293207453,640068.5930263129,639160.7287488431,639384.4084377183,638506.7647683267,637836.7463358919,639303.0205103463,639045.4182576606,637889.5208676661,643608.0395533801,644660.3686917758,644807.5007459201,644055.8020003792,644004.7985236729,643828.7998079469,648002.6180677606,643438.8332405403,643101.3742837426,643068.7117299446,642513.0374975023,643126.1821958068,642312.0254059541,641432.5346495764,644871.4119174794,637187.4700057959,638372.9455688297,643275.946369819,647081.0100076408,642465.5722203177,642889.618991915,640744.1838266278,642744.9053558036,643203.0573592281,640671.3964695348,641060.668531723,639546.3427821912,638954.3595692894,640022.7718857436,639333.0548919217,637951.886678179,639685.7004327765,639272.2176390159,639220.1143686195,639523.9047317485,637811.2030085957,638840.0253393577,644439.6980401423,643708.8236632174,644521.8323119825,645708.5259798458,645662.0286385682]
Upon plotting it looks like this:
[EDIT : after 20min of writing the post to show the results of each plot operation, it turns out that as a new user I can only insert 1 media file, so I will paste them all together at the end of the post]
PyPlot.scatter(xs,ys,s=1)
When using the reprojection function to project to other EPSG (different from 4326) everything works as expected:
transformed_points = reproject_points(xs,ys,2062,3857)
xs_projected = []
ys_projected = []
for i in 1:size(transformed_points)[1]
push!(xs_projected,AG.getx(transformed_points[i],0))
push!(ys_projected,AG.gety(transformed_points[i],0))
end
PyPlot.scatter(xs_projected,ys_projected,s=1)
BUT, and this is the point of this post, when I project to 4326, the x and y coordinates are flipped:
transformed_points = reproject_points(xs,ys,2062,4326)
xs_projected = []
ys_projected = []
for i in 1:size(transformed_points)[1]
push!(xs_projected,AG.getx(transformed_points[i],0))
push!(ys_projected,AG.gety(transformed_points[i],0))
end
PyPlot.scatter(xs_projected,ys_projected,s=1)
Just for convenience, and so it is easy to see that it is the case, let me plot flipping the coordinates, just to show that it plots correctly in that case:
PyPlot.scatter(ys_projected,xs_projected,s=1)
This behaviour is highly undesirable, I do not know whether it emanates from ArchGDAL.jl or from GDAL.jl but it does not allow to make general functions, I would need to start making special cases to deal with EPSG:4326.
Although I know (and cannot understand why) that generally geographers call first lat and then lon, this does not make sense from a software perspective, and I never had this problem in any other implementation on other languages. I think it needs to be corrected, but I have no idea where it comes from and what code in what package would need to change.
Can someone that knows more than me elevate this to the correct package maintainer?
Thank you very much, and if this comes from me not understanding something please let me know.
Carlos.
