# Read and use .gpkg file

**URL:** <https://discourse.julialang.org/t/read-and-use-gpkg-file/117847>\
**Category:** General Usage\
**Tags:** question, geo\
**Created:** [August 5, 2024, 4:30pm UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847 "2024-08-05T16:30:37Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [August 5, 2024, 4:30pm UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/1 "2024-08-05T16:30:37Z")

</div>

Hi,  
Here’s my latest challenge:

I’m looking to get the 5km tiles out of this file: [os\_bng\_grids.7z](https://github.com/OrdnanceSurvey/OS-British-National-Grids/tree/main).

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. 😟

```julia
    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.

```julia
(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!

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 5, 2024, 6:34pm UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/2 "2024-08-05T18:34:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 5, 2024, 7:08pm UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/3 "2024-08-05T19:08:25Z")

</div>

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:

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

```

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

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [August 6, 2024, 7:41am UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/4 "2024-08-06T07:41:18Z")

</div>

> [@juliohm](#):
>
> can you please explain in more detail what you are trying to achieve?

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:

```julia
points = Point.(zip(df.EASTING, df.NORTHING))
newdf = georef(df, points)

```

1. `geojoin` the address points to the grid tile geometry:

```julia
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!)

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [August 6, 2024, 7:51am UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/5 "2024-08-06T07:51:53Z")

</div>

> [@juliohm](#):
>
> 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:
> 
> ```julia
> GeoIO.load("os_bng_grids.gpkg", layer=1)
> 
> ```
> 
> There are 6 layers in this file (0, 1, …, 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

```julia
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.

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [August 6, 2024, 8:07am UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/6 "2024-08-06T08:07:29Z")

</div>

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:

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

```

generates the relatively uninformative message:

```julia
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
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!).

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [August 6, 2024, 1:07pm UTC](https://discourse.julialang.org/t/read-and-use-gpkg-file/117847/7 "2024-08-06T13:07:30Z")

</div>

> [@TimG](#):
>
> 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.

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.

> [@TimG](#):
>
> 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!)

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

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