ESRI code for British National Grid not known by GeoIO

TimG · July 30, 2024, 4:05pm

ERROR: LoadError: ArgumentError: EPSG/ESRI code for the ESRI ID "British_National_Grid" not found in dictionary.
Please check https://github.com/JuliaEarth/CoordRefSystems.jl/blob/main/src/strings.jl
If you know the EPSG/ESRI code of a given ESRI WKT string, please submit a pull request.

I’m reading this shape file from here at the the UK Office for National Statistics (ONS): LAD_DEC_2023_UK_BGC.shp.
Google tells me the British National Grid is EPSG:27700.

How can I circumvent this error and make the file load successfully? Alternatively, can GeoIO be made to “just work” with BNG again as it did a few weeks ago?

nilshg · July 30, 2024, 4:11pm

What is the question?

TimG · July 30, 2024, 4:15pm

Updated OP to clarify.

nilshg · July 30, 2024, 4:22pm

Looking at the code referenced in the error message, it looks like you need to submit a PR that adds "British National Grid" => EPSG{27700} to the esriid2code dictionary.

If your issues are caused by an update of GeoIO, you can either undo that package update if you are still in the same session (see Pkg.undo) or you can downgrade GeoIO to whichever version worked for you by requesting that version specifically, i.e. add GeoIO@x.y where x.y is the version number you want.

juliohm · July 30, 2024, 4:46pm

Thank you @nilshg for the help.

That is precisely it. @TimG you need to find what is the ESRI ID in the CRS string of your file, and submit a PR pointing to the correct EPSG/ESRI code. Please let us know if you need assistance with the PR.

We are gradually updating GeoIO.jl to load CRS information stored in all these file formats. In the past it was just loading assuming a Cartesian CRS.

TimG · July 31, 2024, 7:46am

Thanks both for guidance.

Here is my pull request. I hope it’s right!

Update strings.jl to add British_National_Grid by TimG1964 · Pull Request #132 · JuliaEarth/CoordRefSystems.jl (github.com)

juliohm · July 31, 2024, 9:37am

Thank you! Released a patch with the addition: New version: CoordRefSystems v0.10.3 by JuliaRegistrator · Pull Request #112095 · JuliaRegistries/General · GitHub

The next step consists of mapping the EPSG{27700} code you added to the corresponding type in CoordRefSystems.jl. You can do this by adding a get method here:

github.com

JuliaEarth/CoordRefSystems.jl/blob/67ee7203e52a19b93b7c08017afbe5a15728a6a3/src/get.jl#L29-L46


      
          get(::Type{EPSG{3395}}) = Mercator{WGS84Latest}
          get(::Type{EPSG{3857}}) = WebMercator{WGS84Latest}
          get(::Type{EPSG{4208}}) = LatLon{Aratu}
          get(::Type{EPSG{4618}}) = LatLon{SAD69}
          get(::Type{EPSG{4674}}) = LatLon{SIRGAS2000}
          get(::Type{EPSG{4326}}) = LatLon{WGS84Latest}
          get(::Type{EPSG{4988}}) = Cartesian3D{shift(ITRF{2000}, 2000.4)}
          get(::Type{EPSG{4989}}) = LatLonAlt{shift(ITRF{2000}, 2000.4)}
          get(::Type{EPSG{5527}}) = LatLon{SAD96}
          get(::Type{EPSG{9988}}) = Cartesian3D{ITRF{2020}}
          get(::Type{EPSG{10176}}) = Cartesian3D{IGS20}
          get(::Type{EPSG{32662}}) = PlateCarree{WGS84Latest}
          get(::Type{ESRI{54017}}) = Behrmann{WGS84Latest}
          get(::Type{ESRI{54030}}) = Robinson{WGS84Latest}
          get(::Type{ESRI{54034}}) = Lambert{WGS84Latest}
          get(::Type{ESRI{54042}}) = WinkelTripel{WGS84Latest}
          get(::Type{ESRI{102035}}) = Orthographic{90.0u"°",0.0u"°",true,WGS84Latest}
          get(::Type{ESRI{102037}}) = Orthographic{-90.0u"°",0.0u"°",true,WGS84Latest}

Please let us know if you need assistance.

TimG · July 31, 2024, 9:49am

Yes, please - help!

I have looked here on the EPSG io and there is lots of information but I don’t know what to take from there or how to add it to your table.

Attributes
Unit: metre

Geodetic CRS: OSGB36

Datum: Ordnance Survey of Great Britain 1936

Ellipsoid: Airy 1830

Prime meridian: Greenwich

Data source: EPSG

Information source: Ordnance Survey of Great Britain.

Revision date: 2021-03-18

Scope: Engineering survey, topographic mapping.

Area of use: United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.

Coordinate system: Cartesian 2D CS. Axes: easting, northing (E,N). Orientations: east, north. UoM: m.

juliohm · July 31, 2024, 10:29am

Of course! Let me share some resources here to help you and also help future contributors.

If you know the EPSG code, you can use epsg.org’s search engine to find a more organized set of parameters for your CRS: Geodetic Database

In this case, I searched for 27700 and found your CRS in the CRSs tab:

By clicking on it, you will be directed to the following page:

The first information you need to find is the datum of the CRS. If you start clicking on the links, you will find a datum entry:

In this case the datum is called “Ordnance Survey of Great Britain 1936” or “OSGB36” for short. This is already available in CoordRefSystems.jl as you can see here:

github.com

JuliaEarth/CoordRefSystems.jl/blob/67ee7203e52a19b93b7c08017afbe5a15728a6a3/src/datums.jl#L197-L204


      
          """
              OSGB36
          
          Airy 1830 datum.
          """
          abstract type OSGB36 <: Datum end
          
          ellipsoid(::Type{OSGB36}) = Airy🌎

The next information you need is the coordinate system, which is also there:

It says that you are dealing with Cartesian2D coordinates in m units. This is also available in CoordRefSystems.jl as you can see here:

github.com

JuliaEarth/CoordRefSystems.jl/blob/67ee7203e52a19b93b7c08017afbe5a15728a6a3/src/crs/basic.jl#L53-L70


      
          """
              Cartesian2D(x, y)
              Cartesian2D{Datum}(x, y)
              Cartesian2D((x, y))
              Cartesian2D{Datum}((x, y))
          
          Alias to [`Cartesian`](@ref) with 2 coordinates.
          
          ## Examples
          
          ```julia
          Cartesian2D(1, 1) # add default units
          Cartesian2D(1u"m", 1u"m") # integers are converted converted to floats
          Cartesian2D(1.0u"km", 1.0u"km")
          Cartesian2D{WGS84Latest}(1.0u"m", 1.0u"m")
          ```
          """
          const Cartesian2D{Datum} = Cartesian{Datum,2}

Now you can construct the full CRS type for your EPSG code:

get(::Type{EPSG{27700}}) = Cartesian2D{OSGB36}

Can you submit a PR with this code added, and a test for coverage?

Note

There is a final missing feature in CoordRefSystems.jl to fully handle local grids (e.g. British grid) in conversions. This means that conversions from EPSG{27700} to other CRS may be slightly off. We are planning to add this feature in the following months.

TimG · July 31, 2024, 10:32am

Thank you. I’d never have figured this out alone!

I’m afraid I’ll need a bit more help with this, too.

juliohm · July 31, 2024, 10:41am

Of course, you just need to add a test here besides the other ones:

github.com

JuliaEarth/CoordRefSystems.jl/blob/67ee7203e52a19b93b7c08017afbe5a15728a6a3/test/get.jl#L2-L20


      
          # EPSG/ESRI code
          @test CoordRefSystems.get(EPSG{3395}) === Mercator{WGS84Latest}
          @test CoordRefSystems.get(EPSG{3857}) === WebMercator{WGS84Latest}
          @test CoordRefSystems.get(EPSG{4208}) === LatLon{Aratu}
          @test CoordRefSystems.get(EPSG{4618}) === LatLon{SAD69}
          @test CoordRefSystems.get(EPSG{4674}) === LatLon{SIRGAS2000}
          @test CoordRefSystems.get(EPSG{4326}) === LatLon{WGS84Latest}
          @test CoordRefSystems.get(EPSG{4988}) === Cartesian{CoordRefSystems.shift(ITRF{2000}, 2000.4),3}
          @test CoordRefSystems.get(EPSG{4989}) === LatLonAlt{CoordRefSystems.shift(ITRF{2000}, 2000.4)}
          @test CoordRefSystems.get(EPSG{5527}) === LatLon{SAD96}
          @test CoordRefSystems.get(EPSG{9988}) === Cartesian{ITRF{2020},3}
          @test CoordRefSystems.get(EPSG{10176}) === Cartesian{IGS20,3}
          @test CoordRefSystems.get(EPSG{32662}) === PlateCarree{WGS84Latest}
          @test CoordRefSystems.get(ESRI{54017}) === Behrmann{WGS84Latest}
          @test CoordRefSystems.get(ESRI{54030}) === Robinson{WGS84Latest}
          @test CoordRefSystems.get(ESRI{54034}) === Lambert{WGS84Latest}
          @test CoordRefSystems.get(ESRI{54042}) === WinkelTripel{WGS84Latest}
          @test CoordRefSystems.get(ESRI{102035}) === CoordRefSystems.Orthographic{90.0u"°",0.0u"°",true,WGS84Latest}
          @test CoordRefSystems.get(ESRI{102037}) === CoordRefSystems.Orthographic{-90.0u"°",0.0u"°",true,WGS84Latest}

TimG · July 31, 2024, 10:43am

So just added this:
@test CoordRefSystems.get(::Type{EPSG{227700}}) === Cartesian2D{OSGB36}

TimG · July 31, 2024, 10:50am

I’ve made a new PR but it is against the clone I took earlier this morning which you’ve now changed. Presumably this’ll still work, though.

Updated get.jl and relevant test for EPSG:27700 (British National Grid) by TimG1964 · Pull Request #133 · JuliaEarth/CoordRefSystems.jl (github.com)

juliohm · July 31, 2024, 11:00am

I commented on the PR. We sorted the lines of the dictionary alphabetically and that caused conflicts with your new changes. Do you know how to fix git conflicts? You can click on the “Resolve conflicts” and retain only the newly added lines in get.jl and test / get.jl

juliohm · July 31, 2024, 11:44am

@TimG a small correction regarding the type. In the current design we need to use the type that corresponds to the conversion method:

In this case, this is:

TM = TransverseMercator{0.9996012717,49°,-2°,OSGB36}
get(::EPSG{27700}) = shift(TM, xₒ=400000m, yₒ=-100000m)

Given that this is much more verbose, we will take care of it for you.

Will let you know here when the patch is released.

TimG · July 31, 2024, 11:57am

I just resubmitted my PR as suggested and then came back to read your latest.
Thank you for taking over!

nilshg · July 31, 2024, 12:46pm

Might be useful to put this somewhere in the docs and link to that in the error message rather than straight to the code as you currently do? Granted I’m not a geoscientist, but I don’t think I would have ever figured out to do all these things from the error message…

juliohm · July 31, 2024, 12:55pm

It is in our TODO list. Maybe a link to this discourse thread would already help in the meantime?

juliohm · July 31, 2024, 1:05pm

@TimG new patch released:

github.com/JuliaRegistries/General

New version: CoordRefSystems v0.10.4

JuliaRegistries:master ← JuliaRegistries:registrator-coordrefsystems-b46f11dc-v0.10.4-61402e1cf6

opened 01:04PM - 31 Jul 24 UTC

JuliaRegistrator

+3 -0

- Registering package: CoordRefSystems - Repository: https://github.com/JuliaEar…th/CoordRefSystems.jl - Created by: @juliohm - Version: v0.10.4 - Commit: 1eb518ca50881db8203c660f77ec045f868a1ed3 - Reviewed by: @juliohm - Reference: https://github.com/JuliaEarth/CoordRefSystems.jl/commit/1eb518ca50881db8203c660f77ec045f868a1ed3#commitcomment-144854404 - Description: Unitful coordinate reference systems for geographic maps in Julia

Please let us know if you are still unable to load the shapefile with correct CRS.

TimG · July 31, 2024, 1:25pm

Thank you for turning this round so quickly.

However:

ERROR: LoadError: MethodError: no method matching topoint(::Shapefile.Point, ::Type{CoordRefSystems.ShiftedCRS{CoordRefSystems.TransverseMercator{0.9996012717, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(49.0), Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(-2.0), OSGB36}, Quantity{Float64, Unitful.Dimensions{()}(), Unitful.FreeUnits{(Unitful.Unit{:Degree, Unitful.Dimensions{()}()}(0, 1//1),), Unitful.Dimensions{()}(), nothing}}(0.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(400000.0), Quantity{Float64, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1),), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}(), nothing}}(-100000.0)}})

Closest candidates are:
  topoint(::Any, ::Type{<:GeodeticLatLon{Datum}}) where Datum
   @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:62
  topoint(::Any, ::Type{<:Cartesian3D{Datum}}) where Datum
   @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:60
  topoint(::Any, ::Type{<:Cartesian2D{Datum}}) where Datum
   @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:58

Stacktrace:
  [1] (::GeoIO.var"#8#9"{UnionAll})(p::Shapefile.Point)
    @ GeoIO .\none:0
  [2] iterate
    @ .\generator.jl:47 [inlined]
  [3] collect(itr::Base.Generator{Base.Generator{UnitRange{Int64}, GeoInterface.var"#22#23"{GeoInterface.LinearRingTrait, Shapefile.LinearRing{Shapefile.Point, Nothing, Nothing}}}, GeoIO.var"#8#9"{UnionAll}})
    @ Base .\array.jl:834
  [4] topoints(geom::Shapefile.LinearRing{Shapefile.Point, Nothing, Nothing}, CRS::Type)
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:64
  [5] tochain(geom::Shapefile.LinearRing{Shapefile.Point, Nothing, Nothing}, CRS::Type)
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:67
  [6] (::GeoIO.var"#toring#10"{UnionAll})(g::Shapefile.LinearRing{Shapefile.Point, Nothing, Nothing})
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:81
  [7] topolygon(geom::Shapefile.SubPolygon{Shapefile.LinearRing{Shapefile.Point, Nothing, Nothing}}, CRS::Type)
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:82
  [8] #13
    @ .\essentials.jl:0 [inlined]
  [9] iterate
    @ .\generator.jl:47 [inlined]
 [10] collect(itr::Base.Generator{Vector{Shapefile.SubPolygon{Shapefile.LinearRing{Shapefile.Point, Nothing, Nothing}}}, GeoIO.var"#13#14"{UnionAll}})
    @ Base .\array.jl:834
 [11] togeometry(::GeoInterface.MultiPolygonTrait, geom::Shapefile.Polygon, crs::GeoFormatTypes.ESRIWellKnownText{GeoFormatTypes.CRS})
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:108
 [12] geom2meshes
    @ C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:112 [inlined]
 [13] geom2meshes
    @ C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\conversion.jl:111 [inlined]
 [14] _broadcast_getindex_evalf
    @ .\broadcast.jl:709 [inlined]
 [15] _broadcast_getindex
    @ .\broadcast.jl:682 [inlined]
 [16] getindex
    @ .\broadcast.jl:636 [inlined]
 [17] copy
    @ .\broadcast.jl:942 [inlined]
 [18] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(GeoIO.geom2meshes), Tuple{Vector{Union{Missing, Shapefile.Polygon}}, Base.RefValue{GeoFormatTypes.ESRIWellKnownText{GeoFormatTypes.CRS}}}})
    @ Base.Broadcast .\broadcast.jl:903
 [19] asgeotable(table::Shapefile.Table{Union{Missing, Shapefile.Polygon}})
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\utils.jl:13
 [20] load(fname::String; repair::Bool, layer::Int64, kwargs::@Kwargs{})
    @ GeoIO C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\load.jl:101
 [21] load
    @ C:\Users\TGebbels\.julia\packages\GeoIO\V1sFL\src\load.jl:26 [inlined]
 [22] macro expansion
    @ .\timing.jl:279 [inlined]
 [23] main()
    @ Main c:\Users\TGebbels\OneDrive - The National Lottery Community Fund\Documents\DCMS Database\CityMaps\GeoStat City Maps.jl:166
 [24] top-level scope
    @ c:\Users\TGebbels\OneDrive - The National Lottery Community Fund\Documents\DCMS Database\CityMaps\GeoStat City Maps.jl:206
 [25] include(fname::String)
    @ Base.MainInclude .\client.jl:489
 [26] run(debug_session::VSCodeDebugger.DebugAdapter.DebugSession, error_handler::VSCodeDebugger.var"#3#4"{String})
    @ VSCodeDebugger.DebugAdapter c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\DebugAdapter\src\packagedef.jl:126
 [27] startdebugger()
    @ VSCodeDebugger c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.105.2\scripts\packages\VSCodeDebugger\src\VSCodeDebugger.jl:45
 [28] top-level scope
    @ c:\Users\TGebbels\.vscode\extensions\julialang.language-julia-1.105.2\scripts\debugger\run_debugger.jl:12
 [29] include(mod::Module, _path::String)
    @ Base .\Base.jl:495
 [30] exec_options(opts::Base.JLOptions)
    @ Base .\client.jl:318
 [31] _start()
    @ Base .\client.jl:552
in expression starting at c:\Users\TGebbels\...\Documents\DCMS Database\CityMaps\GeoStat City Maps.jl:206