Converting Longitude & Latitude to X Y on a map using Julia


#1

Hey. I searched for this question and I found this link: https://stackoverflow.com/questions/1369512/converting-longitude-latitude-to-x-y-on-a-map-with-calibration-points
My question is not exactly the same as the above link, because I do not have any calibration points at all. All I have is a bunch of locations in the latitudes and longitudes, and I need to project them to a self-designed planar map. In fact, I think if the concerned area is not big enough, then maybe I can use the latitude and longtitude directly as the x axis and y axis of the points in the planar map.
Any ideas or comments about this simple task?


#2

There are other options (like using proj4.jl) but the one I know is with GMT.jl. With its mapproject module you would do it like this (the [0 0;10 10;20 20;30 30] is an invented lon lat array). Here I’m projecting from geographic to mercator. The syntax after the -J is that of the proj4 lib (ah, and need GMT6dev here, otherwise the GMT projection syntax must be used).

julia> gmt("mapproject -J+proj=merc", [0 0;10 10;20 20;30 30])
1-element Array{GMT.GMTdataset,1}:
 GMT.GMTdataset([0.0 7.08115e-10; 1.11319e6 1.11148e6; 2.22639e6 2.25842e6; 3.33958e6 3.48219e6], Any[], "", Any[], "", "")

#3

So, you are projecting four locations with lon lat pair (0,0), (10,10), (20,20), (30,30) to four points (0.0, 7.08115e-10), (1.11319e6, 1.11148e6), (2.22639e6, 2.25842e6), (3.33958e6, 3.48219e6)?


#4

@joa-quim, is there any advantage to do this? What are the map size and scale?


#5

I don’t know exactly what you want/need to do. My example just converted a couple of invented points to Mercator meters. There is no scale it nor map size. The converted xx mean meters east of Greenwich and yy meters north of equator. For a more realistic case maybe an UTM projection would be better, but as I said I don’t know what for you need the projected coordinates.


#6

If the area is not too large (thousands of kilometers), a simple way is to convert to UTM coordinates within the same zone. UTM consists of an x coordinate (easting) and y coordinate (northing), and a change in 1 in easting or northing is roughly equivalent to a distance of 1 meter. This makes the coordinates conceptually easy to understand, and it becomes trivial to calculate approximate distances. You can use the Geodesy package for this. An example:

using Geodesy

# latitude and longitude for two cities in Sweden
stockholm_lla = LLA(59.3293, 18.0686)
uppsala_lla = LLA(59.8586, 17.6389)

# all of Sweden uses UTM zone 33 north
utm_sweden = UTMfromLLA(33, true, wgs84)

Now to calculate the UTM coordinates:

julia> stockholm_utm = utm_sweden(stockholm_lla)
UTM(674571.8664470176, 6.580743008454201e6, 0.0)

julia> uppsala_utm = utm_sweden(uppsala_lla)
UTM(647793.4861253307, 6.638608056491923e6, 0.0)

julia> diff = stockholm_utm.x - uppsala_utm.x, stockholm_utm.y - uppsala_utm.y
(26778.38032168697, -57865.048037721775)

So Uppsala is roughly 27 km west and 58 km north of Stockholm, for a total distance of

julia> √sum(diff.^2)
63760.84564261018

approximately 64 km. Let’s compare that with a more accurate distance calculation to see how bad the approximation is:

julia> distance(stockholm_utm, uppsala_utm, 33, true)
63765.736076373265

julia> distance(stockholm_lla, uppsala_lla)
63765.73607637202

Off by ~5 meters, which is good enough for many applications (and several orders of magnitude faster to calculate).


#7

Very good. Thank you so much. Was just curious for the following two questions.
(1): Why all of Sweden use UTM zone 33 north, does this have something to do with the time zone of Sweden?
(2): The last two distance function calls, what are the algorithm do they use to calculate the more accurate distance?


#8

One more question which is more revenant. If I plan to visualise it in a map, it’s probably good to have the size of the map, something like xMin, xMax, yMin, yMax. So, I think the boundary of all my points would roughly correspond to the approximated boundary of the map.


#9

No, time zone has nothing to do with this. Geographically, Sweden extend over zones 32 - 35, but the main part is in zone 33, so the national projection uses that zone. It leads to some distorsions towards the edges, but it’s much easier to have a single zone.

To figure that out, you can either read the docs, or use the @less macro to check the source:

julia> @less distance(stockholm_lla, uppsala_lla)
distance(a::LLA, b, datum = wgs84) = distance(ECEFfromLLA(datum)(a), b, datum)
...
julia> @less distance(ECEFfromLLA(wgs84)(stockholm_lla), ECEFfromLLA(wgs84)(uppsala_lla))
distance(a::T, b::T) where {T <: AbstractVector} = norm(a-b)

So it’s cartesian distance in Earth-Centered-Earth-Fixed (ECEF) coordinates (not great circle distance).

What’s the question here? :slight_smile: