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


Hey. I searched for this question and I found this link:
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?


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[], "", "")


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)?


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


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.


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)

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)

julia> distance(stockholm_lla, uppsala_lla)

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


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?


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.


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: