My Blog - Geodetic coordinate transforms with Julia

This blog post is about geodetic coordinate transforms using the Julia programming language.

In applied geophysics, the following question comes up frequently:

Problem 💻

How can I transform geographical coordinates given as “latitude, longitude” into some other coordinate system, e.g., UTM?

The Julia package Proj.jl offers all functionality that is required.

After installing the package using Julia’s package manager from the REPL, we are ready to go:

using Proj
using DataFrames, CSV

Let’s assume that we have downloaded a set of coordinates from a handheld GPS receiver. The content of the data file coords.dat may look like this:

# lat, lon
50.924833, 13.330562
50.982648, 13.530406

We first read in the data:

df = DataFrame(CSV.File("coords.dat", delim=","))

2 rows × 2 columns

	# lat	lon
	Float64	Float64
1	50.9248	13.3306
2	50.9826	13.5304

Next we arrange the data such that it is suitable for processing with Proj.jl:

latlon = Array(df)

2×2 Matrix{Float64}:
 50.9248  13.3306
 50.9826  13.5304

The following step is essential. Since we transform data from one coordinate system into another, we have to inform Proj.jl about the source and target systems. To this end, we exploit the convenient EPSG codes.

trans = Proj.Transformation("EPSG:4326", "EPSG:25833")

Transformation
    source: WGS 84
    target: ETRS89 / UTM zone 33N
    direction: forward

The next lines will finally transform our GPS coordinates into UTM zone 33 coordinates which we refer to as easting and northing:

etrs = [trans(latlon[i, :]) for i in 1:size(latlon, 1)]
UTM = hcat(collect.(etrs)...);

easting = UTM[1, :]
northing = UTM[2, :]
@show easting
@show northing;

easting = [382670.42457542894, 396843.0546304857]
northing = [5.642793241297958e6, 5.648923257065746e6]