# So now I convert my OSGB36 data to LLA ... oh ... that looks ... hard

Hi,

The UK’s Office of National Statistics publishes data which uses OSGB36 East / North coordinates.

Of which I have 39 million addresses from the ONSUD database

e.g. East/North “426642,380231”

Using an online converter

this is

`Lat:53.3183427137348 Long:-1.60153684101643`

The other data I have is all in LL

I don’t care too much about the direction of conversion so long as I can use a common format

I tried using Geodesy as this seems to be the tool to use but I don’t understand the results

``````
julia> ecef = ECEFfromLLA(osgb36)(LLA(53.3183427137348, -1.60153684101643))
ECEF(3.8164622703267145e6, -106705.89246363874, 5.09141104371438e6)

julia> @printf("(%d, %d, %d)", ecef...)
(3816462, -106706, 5091411)
``````

But that makes no sense, clearly I have much to learn.

You can use `Proj4.jl` to convert your OSGB36 X, Y data to WGS84 Lat, Lon:

``````using Proj4
src = Projection(Proj4.epsg[27700])     #  [m] OSGB 1936 / British National Grid
tgt = Projection(Proj4.epsg[4326])      # [deg] WGS84 = ETRS89 = EPSG: 4326
Xs, Ys = 426642.0, 380231.0
Proj4.transform(src, tgt, [Xs, Ys])     #  -1.6015368, 53.3183427
``````
2 Likes

You can also do this with `reproject` in ArchGDAL.jl

1 Like

thank you both

1 Like

If I understand correctly, with `ECEFfromLLA(osgb36)(LLA(x, y))`, what you are doing is to convert latitude longitude coordinate with no elevation, that lie on the OSGB36 datum, to a Earth-centered, Earth-fixed (ECEF) coordinate system, i.e. measured from the center of the Earth. What you want to do instead is to convert coordinates that you have that are in the OSGB 1936 / British National Grid (EPSG:27700) coordinate system, to latitudes and longitudes.

As mentioned,

However I’d strongly recommend to use the new API, both for increased accuracy, and because the old API is about to be removed:

``````julia> using Proj4

julia> trans = Proj4.Transformation("EPSG:27700", "EPSG:4326", always_xy=true)
Transformation
source: OSGB 1936 / British National Grid
target: WGS 84 (with axis order normalized for visualization)

julia> trans((426642.0, 380231.0))
2-element SVector{2, Float64} with indices SOneTo(2):
-1.6015358313331525
53.3183599763518
``````

This removal will happen as part of the renaming to Proj.jl in Rename package to Proj.jl by visr · Pull Request #57 · JuliaGeo/Proj4.jl · GitHub. Besides renaming the package, this example will continue to work.

1 Like

Ah yes, the no-elevation explanation makes sense.

Curiously I get a different result to your example

``````julia> trans = Proj4.Transformation("EPSG:27700", "EPSG:4326", always_xy=true)
Transformation
source: OSGB 1936 / British National Grid
target: WGS 84 (with axis order normalized for visualization)

julia> trans((426642.0, 380231.0))
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
-1.601536798482117
53.31834270000983
``````

And both are also different to the previous solution

``````julia> Proj4.transform(Projection(Proj4.epsg[27700]), Projection(Proj4.epsg[4326]), [426642.0, 380231.0])
2-element Vector{Float64}:
-1.601536844755606
53.31834269853471
``````

but for my usage that is far enough down the decimal to not make too much difference, as it is in the 1m realm.

Hi,

an FYI for anyone reading later

this new API is not threadsafe

I replaced this

`en2lalo(e, n) = Proj4.transform(Projection(Proj4.epsg[27700]), Projection(Proj4.epsg[4326]), [e, n])`
with this
`en2lalo(e, n) = Proj4.Transformation("EPSG:27700", "EPSG:4326", always_xy=true)((e,n))`

And using mutliple threads (eleven) I get a `received signal: 11` error

Commenting out the `Threads.@threads` made it work again

``````##Threads.@threads
end
``````

I’ve raised an Issue with the project - Thread Safety - Signal 11 faults if using multiple threads · Issue #58 · JuliaGeo/Proj4.jl · GitHub

1 Like

Accuracy in the sixth decimal place in the latitude measurement is worth about 10 cm.
It is almost certain that common input data is far from that kind of accuracy.

4 Likes

creating one en2lalo function per thread solved this issue

I was scratching my head why I had slightly different results with GMT (that uses GDAL for this calculation). It turns out that is not indifferent if the EPSG code is converted into aPROJ4 string or a WKT. (I’m using the `println` to force writing with higher number of decimals)

``````println(xy2lonlat([426642.0 380231.0], s_srs=epsg2wkt(27700)))
[-1.601536798482117 53.31834270000983]

println(xy2lonlat([426642.0 380231.0], s_srs=epsg2proj(27700)))
[-1.6000274390369076 53.31806612704705]
``````

former result is the same as we get when using the CLI `cs2cs`

``````echo 426642.0 380231.0 | cs2cs EPSG:27700 EPSG:4326 -f %.12f
53.318342700010 -1.601536798482 0.000000000000
``````

so it must be the most accurate one.
It’s the first time I see such a sensible difference between the use of `PROJ4` and `WKT` strings.