Read and use .gpkg file

General Usage
,
1

Hi,
Here’s my latest challenge:

I’m looking to get the 5km tiles out of this file: os_bng_grids.7z.

I extracted the content, a 200MB gpgk file of the same name as the 7z file, and can load this effortlessly using GeoIO.jl but I cannot coax it to reveal its structure or content in a way that I can access - or at all infact. :worried:

    geo = GeoIO.load("os_bng_grids.gpkg")
    println(values(geo))
    println((geo[geo.tile_name.=="HL", ["geometry"]]))

shows the tile_names, but these are the names of the 100km squares. I’m looking for the names and grid coords of the 5km squares.

(tile_name = ["HL", "HM", "HN", "HO", "HP", "HQ", "HR", "HS", "HT", "HU", "HV", "HW", "HX", "HY", "HZ", "JL", "JM", "JQ", "JR", "JV", "JW", "NA", "NB", "NC", "ND", "NE", "NF", "NG", "NH", "NJ", "NK", "NL", "NM", "NN", "NO", "NP", "NQ", "NR", "NS", "NT", "NU", "NV", "NW", "NX", "NY", "NZ", "OA", "OB", "OF", "OG", "OL", "OM", "OQ", "OR", "OV", "OW", "SA", "SB", "SC", "SD", "SE", "SF", "SG", "SH", "SJ", "SK", "SL", "SM", "SN", "SO", "SP", "SQ", "SR", "SS", "ST", "SU", "SV", "SW", "SX", "SY", "SZ", "TA", "TB", "TF", "TG", "TL", "TM", "TQ", "TR", "TV", "TW"],)
1×1 GeoTable over 1 view(::GeometrySet, [1])

This must be simple, but it’s currently defeating me!

Can anyone point me in the right direction?

Thanks!

2

Hi @TimG , can you please explain in more detail what you are trying to achieve? Please tag your questions with geo or submit it to the Geo community here on Discourse. We don’t get notifications otherwise.

3

I think I got you. What you need is to specify the layer of the file that you want to load. The package loads layer=0 by default. Try the following:

GeoIO.load("os_bng_grids.gpkg", layer=1)

There are 6 layers in this file (0, 1, …, 5).

1 Like
4

I have a largish dataset of UK addresses in a dataframe. Each address has a grid reference locating it to the nearest metre on the British National Grid. I need to tag each address with the name of the 5km grid tile it lies in.

If/when I’ve read and understood the grid data, my plan of attack was

  1. Convert the grid reference to a geometry, like:
points = Point.(zip(df.EASTING, df.NORTHING))
newdf = georef(df, points)
  1. geojoin the address points to the grid tile geometry:
joineddf = geojoin(newdf, tilegeometry, kind=:left, pred=((g1, g2) -> intersects(g1, g2)))

This may be a bit of a sledge hammer to crack a nut, though. The grid tiles are systematically laid out. The key bit I don’t know is the two letter codes used for the 100km tiles (which are also the first two characters of the 5km tile names). Once the two letters are known, the digits are easily computable from the full grid reference, I think.

So I am now thinking that I just need to find the co-ordinates of the south-west corner of each 100km grid tile and the associated two-letter code for that tile. From that I can create a look-up from any grid location to the relevant 5km grid tile. I only need to do this once and then using the look-up ought to be faster than using geojoin (although geojoin seems quite fast and my code is generally very inefficient!)

1 Like
5

Thank you. Yes, the layers are the key that I didn’t know. Layer 4 contains the 5km grid and has the following column names:
["tile_name", "1km_grid_ref", "geometry"]

The documentation for this dataset says

The 5km grid also contains '1km_grid_ref' which relates to the 1km grid cell
at the South-West corner of each cell
 e.g. 'tile_name' = "SP19SW" and 'tile_grid_ref' = "SP1090"

I actually need the name of the tile_grid_ref, which is made up of the two-letter code for the 100km tile it lies within and the 4-digit grid reference of the south-west corner of the relevant 5km tile within it.

6

In what is presumably a bad habit, I don’t develop my code in the REPL, but use VSCode, and just run without debugging as I go along.

Doing this to load the gpkg file:

geo = GeoIO.load("os_bng_grids.gpkg", layer=4)
println(geo)
println((geo[1:10, ["geometry"]]))

generates the relatively uninformative message:

36400×3 GeoTable over 36400 GeometrySet
10×1 GeoTable over 10 view(::GeometrySet, 1:10)

which is a little too much of a black box for me!

If I do the same in the REPL, I get

julia> geo = GeoIO.load("os_bng_grids.gpkg", layer=4)
36400×3 GeoTable over 36400 GeometrySet
┌─────────────┬──────────────┬────────────────────────────────────────────────────────────────────────┐
│  tile_name  │ 1km_grid_ref │                                geometry                                │
│ Categorical │ Categorical  │                                PolyArea                                │
│  [NoUnits]  │  [NoUnits]   │                  🖈 ShiftedTransverseMercator{OSGB36}                    │
├─────────────┼──────────────┼────────────────────────────────────────────────────────────────────────┤
│   HL00NE    │    HL0505    │ PolyArea((x: 5000.0 m, y: 1.205e6 m), ..., (x: 5000.0 m, y: 1.21e6 m)) │
│   HL00NW    │    HL0005    │    PolyArea((x: 0.0 m, y: 1.205e6 m), ..., (x: 0.0 m, y: 1.21e6 m))    │
│   HL00SE    │    HL0500    │ PolyArea((x: 5000.0 m, y: 1.2e6 m), ..., (x: 5000.0 m, y: 1.205e6 m))  │
│   HL00SW    │    HL0000    │    PolyArea((x: 0.0 m, y: 1.2e6 m), ..., (x: 0.0 m, y: 1.205e6 m))     │
│   HL01NE    │    HL0515    │ PolyArea((x: 5000.0 m, y: 1.215e6 m), ..., (x: 5000.0 m, y: 1.22e6 m)) │
│   HL01NW    │    HL0015    │    PolyArea((x: 0.0 m, y: 1.215e6 m), ..., (x: 0.0 m, y: 1.22e6 m))    │
│   HL01SE    │    HL0510    │ PolyArea((x: 5000.0 m, y: 1.21e6 m), ..., (x: 5000.0 m, y: 1.215e6 m)) │
│   HL01SW    │    HL0010    │    PolyArea((x: 0.0 m, y: 1.21e6 m), ..., (x: 0.0 m, y: 1.215e6 m))    │
│   HL02NE    │    HL0525    │ PolyArea((x: 5000.0 m, y: 1.225e6 m), ..., (x: 5000.0 m, y: 1.23e6 m)) │
│   HL02NW    │    HL0025    │    PolyArea((x: 0.0 m, y: 1.225e6 m), ..., (x: 0.0 m, y: 1.23e6 m))    │
│   HL02SE    │    HL0520    │ PolyArea((x: 5000.0 m, y: 1.22e6 m), ..., (x: 5000.0 m, y: 1.225e6 m)) │
│   HL02SW    │    HL0020    │    PolyArea((x: 0.0 m, y: 1.22e6 m), ..., (x: 0.0 m, y: 1.225e6 m))    │
│   HL03NE    │    HL0535    │ PolyArea((x: 5000.0 m, y: 1.235e6 m), ..., (x: 5000.0 m, y: 1.24e6 m)) │
│   HL03NW    │    HL0035    │    PolyArea((x: 0.0 m, y: 1.235e6 m), ..., (x: 0.0 m, y: 1.24e6 m))    │
│   HL03SE    │    HL0530    │ PolyArea((x: 5000.0 m, y: 1.23e6 m), ..., (x: 5000.0 m, y: 1.235e6 m)) │
│   HL03SW    │    HL0030    │    PolyArea((x: 0.0 m, y: 1.23e6 m), ..., (x: 0.0 m, y: 1.235e6 m))    │
│   HL04NE    │    HL0545    │ PolyArea((x: 5000.0 m, y: 1.245e6 m), ..., (x: 5000.0 m, y: 1.25e6 m)) │
│   HL04NW    │    HL0045    │    PolyArea((x: 0.0 m, y: 1.245e6 m), ..., (x: 0.0 m, y: 1.25e6 m))    │
│   HL04SE    │    HL0540    │ PolyArea((x: 5000.0 m, y: 1.24e6 m), ..., (x: 5000.0 m, y: 1.245e6 m)) │
│   HL04SW    │    HL0040    │    PolyArea((x: 0.0 m, y: 1.24e6 m), ..., (x: 0.0 m, y: 1.245e6 m))    │
│   HL05NE    │    HL0555    │ PolyArea((x: 5000.0 m, y: 1.255e6 m), ..., (x: 5000.0 m, y: 1.26e6 m)) │
│      ⋮      │      ⋮       │                                   ⋮                                    │
└─────────────┴──────────────┴────────────────────────────────────────────────────────────────────────┘
                                                                                     36379 rows omitted

julia>

which is much more informative! (Wouldn’t have told me about the layers, but does open the black-box dataset a bit!).

1 Like
7

Are you aware of the VSCode shortcuts to run lines of the script in the REPL? For example, press Ctrl+Enter to send the line to the REPL. You can press Ctrl+Shift+p and type “Julia” to learn about all commands related to Julia.

You are gradually acquiring new geospatial skills, and that is really nice! :tada:

These functions will get faster over time. There are various performance optimizations in our TODO list related to geojoin.

1 Like