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

Specific Domains Geo
, ,
1

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.

2

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
3

You can also do this with reproject in ArchGDAL.jl

1 Like
4

thank you both

1 Like
5

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/Proj.jl · GitHub. Besides renaming the package, this example will continue to work.

1 Like
6

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.

7

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
    for i in 1:length(readers)
        dbs[i] = generate(readers[i])
    end

I’ve raised an Issue with the project - https://github.com/JuliaGeo/Proj4.jl/issues/58

1 Like
8

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.

9
4 Likes
10

creating one en2lalo function per thread solved this issue

11

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.