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

General Usage
, ,
1

Hey. I searched for this question and I found this link: gps - Converting Longitude & Latitude to X Y on a map with Calibration points - Stack Overflow
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?

1 Like
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[], "", "")
2 Likes
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).

6 Likes
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?

1 Like
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:

2 Likes
10

Hello,

I was looking for an answer to a similar question and found python - Converting x, y, zoom for Calculating tile bounds - Geographic Information Systems Stack Exchange where it was suggested to read up on web mercator projection. I’m not a specialist of this field but I’m wondering if there is any (pure) Julia library which implements this.
Does Geodesy.jl provides that ? Ping @garborg @c42f

Kind regards

11

I don’t know about a Julia library, but webmercator projections are straightforward. This code is in a non-Julia language, but should be obvious to translate. Latitude and longitude in WGS84

// WGS84 earth semi-minor axis in meters.
R = 6378137.0
easting, northing = R*longRads, R*math.Log(math.Tan(math.Pi*0.25+0.5*latRads))

longRads = easting/R
latRads = 2.0*(math.Atan(math.Exp(northing/R))  - math.Pi*0.25)

Reference:
http://earth-info.nga.mil/GandG/wgs84/web_mercator/(U)%20NGA_SIG_0011_1.0.0_WEBMERC.pdf
equations (4) & (5.a), p. A-2.

1 Like
12

Not a pure Julia library but just in case: Proj4.jl wraps the PROJ cartographic projections library, which includes Web Mercator EPSG:3857 among many others.

1 Like
13

I think a nice projection package would be great. These days I would probably just use LibPQ.jl for a PostGIS database and make most of the work in the PostgreSQL server and get it back to Julia.

2 Likes
14

Which, under the hood, would be using the same PROJ4 external lib that Proj.jl wraps.

3 Likes
15

Aye. I like the database component especially within a larger project (e.g., GUI like pgAdmin/Compass are nice for quick interactive rendering of maps and such). For simple projections and transformations I think Proj4.jl does the job really conveniently.

16

The PlotShapefiles package has a number of functions to help drawing in web mercator projection and overlaying plots onto a map:

See here for these coordinate conversions:

  • WGS84 longitude & latitude to Web Mercator (which is used by Google maps and most web map apps)
  • Web Mercator back to longitude and latitude
  • Longitude, latitude and zoom level to pixel numbers (specifically for the Google API)

And here are some helper functions to overlay onto a Google map
E.g.

  • Given a bounding box and zoom, generating the right Google API info (centre coords, width and height)
  • Building the Google API string
2 Likes
17

Thanks @Wedg . Your project is interesting.

I had a look at your code… maybe you should define some constants and also use some functions (such as deg2rad) to avoid to repeat code, constants.

const POLAR_RADIUS = 6378137.0
const MAX_LAT = rad2deg(2*atan(ℯ^pi)-pi/2)  # 85.051129...
18

Thanks Sébastien (@scelles )
Yes agreed - that would help readability and avoiding repetition. Do you know if there is any performance loss to defining constants outside function definitions (vs defining them inside the function)?

19

If you define them (globally) as const I think it doesn’t hurt. But I’m not a specialist of Julia code optimization… but I ever noticed a performance loss with global variable which wasn’t defined as const.

20

Declare them const there will be no performance hit. If you forget this there will be a big performance cost because the compiler won’t be able to know the type or value of the global variable and will need to generate dynamic non-optimized code for every calculation involving the global variable.

Just a small note @scelles - if you could start a new discussion for new questions that would be better than tacking it onto an old thread like this one.

Regarding your question, Web Mercator is very simple so you could just implement it as already suggested in this thread. TBH we should probably have Web Mercator support in Geodesy.jl because it’s ubiquitous (due to its extreme simplicity — in other mathematical technical aspects it’s quite bad :laughing: )

3 Likes