# Accessing x, y values in a geopackage

**URL:** <https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196>\
**Category:** Geo\
**Created:** [August 25, 2023, 4:45pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196 "2023-08-25T16:45:31Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![chadagreene](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chadagreene/32/34924_2.png) [@chadagreene](https://discourse.julialang.org/u/chadagreene)\
**Post date:** [August 25, 2023, 4:45pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/1 "2023-08-25T16:45:31Z")

</div>

Basic question: I can load a geopackage with:

```julia
using GeoTables 

table = GeoTables.load("myfile.gpkg")

```

Now how do I access the x, y values in the table?

---

<div class="post-metadata">

**Author:** ![alex-s-gardner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex-s-gardner/32/30210_2.png) [@alex-s-gardner](https://discourse.julialang.org/u/alex-s-gardner)\
**Post date:** [August 25, 2023, 4:52pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/2 "2023-08-25T16:52:31Z")

</div>

Would you mind posting a link to an example geopackage?

---

<div class="post-metadata">

**Author:** ![chadagreene](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chadagreene/32/34924_2.png) [@chadagreene](https://discourse.julialang.org/u/chadagreene)\
**Post date:** [August 25, 2023, 4:53pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/3 "2023-08-25T16:53:24Z")

</div>

@alex-s-gardner [https://download.geoservice.dlr.de/icelines/files/Ross1/Ross1\_monthly.zip](https://download.geoservice.dlr.de/icelines/files/Ross1/Ross1_monthly.zip)

---

<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 25, 2023, 5:01pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/4 "2023-08-25T17:01:57Z")

</div>

You can access the column of geometries as if it was a DataFrame @chadagreene:

```julia
geotable.geometry

```

The entire DataFrame syntax is supported:

[https://juliaearth.github.io/GeoStats.jl/stable/quickstart.html](https://juliaearth.github.io/GeoStats.jl/stable/quickstart.html)

Worth checking the quickstart above.

---

<div class="post-metadata">

**Author:** ![chadagreene](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chadagreene/32/34924_2.png) [@chadagreene](https://discourse.julialang.org/u/chadagreene)\
**Post date:** [August 25, 2023, 5:16pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/5 "2023-08-25T17:16:12Z")

</div>

Thanks @juliohm, but it still feels like there’s a barrier between me and the coordinates in the table. How would you complete these lines?

```julia
using GeoTables 

table = GeoTables.load("myfile.gpkg")
x = ???
y = ???

```

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [August 25, 2023, 5:32pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/6 "2023-08-25T17:32:46Z")

</div>

Try

using GMT

D = gmtread(“file.gpkg”)

And the coordinates (and the rest) will be plainly visible

---

<div class="post-metadata">

**Author:** ![alex-s-gardner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex-s-gardner/32/30210_2.png) [@alex-s-gardner](https://discourse.julialang.org/u/alex-s-gardner)\
**Post date:** [August 25, 2023, 5:50pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/7 "2023-08-25T17:50:28Z")

</div>

@juliohm I took a quick crack at this and it really isn’t intuitive how to do this even as someone that lives in Julia, let alone a newcomer… here’s my thought process:

```julia
using DataFrames
using GeoTables
using Downloads

Downloads.download("https://download.geoservice.dlr.de/icelines/files/Ross1/Ross1_monthly.zip", "Ross1_monthly.zip");
run(`unzip $(pwd())/Ross1_monthly.zip`)

```

```julia
table = GeoTables.load("fronts/1SDH_201706-Ross1.gpkg");
# let look at the type
println(typeof(table));
Meshes.MeshData{Meshes.GeometrySet{2,Float64,Meshes.MultiRope{2,Float64,Meshes.Rope{2,Float64,Vector{Meshes.Point2}}}},Dict{Int64,NamedTuple{(:DATE_, :name, :updated, :version, :s1name),NTuple{5,Vector{String}}}}}

```

oh man that looks confusing

but if you look at the display in the REPL you get useful info

```julia
table
1 GeometrySet{2,Float64}
  variables (rank 1)
    └─DATE_ (String)
    └─name (String)
    └─updated (String)
    └─version (String)
    └─s1name (String)

```

but I don’t see geomtery  
lets check the names in table

```julia
names(table)
6-element Vector{String}:
 "DATE_"
 "name"
 "updated"
 "version"
 "s1name"
 "geometry"

```

There it is

```julia
table.geometry
1 GeometrySet{2,Float64}
  └─MultiRope{2,Float64}

```

let’s index into table.geometry

```julia
table.geometry[1]
MultiRope{2,Float64}
  └─Rope(Point(-351057.6779158652, -1.1820828635592712e6), Point(-351066.50089631113, -1.1820828635592712e6), Point(-351066.50089631113, -1.1822428635592712e6)....

```

OK so it’s made up of a bunch of points, I assume these are X,Y  
next I’ll try to index into table.geometry[1] to get a single point

```julia
table.geometry[1][1]
ERROR: MethodError: no method matching getindex(::Meshes.MultiRope{2, Float64, Meshes.Rope{2, Float64, Vector{Meshes.Point2}}}, ::Int64)
Stacktrace:
 [1] top-level scope
   @ REPL[33]:1

```

hmm… I’m a bit stuck now… I’ll go to the GeoTable documentation to see if I can find a function to retrieve point locations

---

<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 25, 2023, 6:19pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/8 "2023-08-25T18:19:45Z")

</div>

Thank you @alex-s-gardner for sharing the experience. We are improving it soon with more helpful displays. Let me go over the details quickly here. I am in the middle of a short-course so can’t type much…

The geometry column has geometries, if you need to access information about the geometries you can use the `@transform`:

```julia
newgeotable = @transform(geotable, :center = centroid(:geometry))

```

If you need to access specific coordinates you do it again:

```julia
@transform(newgeotable, :X = first(coordinates(:center))

```

You can construct as many columns as needed with any Julia expression using the `@transform` macro. Worth reading the split-apply-combine section of the docs.

These macros will preserve the geospatial information all the way, so you can always send the result to the `viewer` at any time:

```julia
geotable |> viewer

```

You can also avoid the macros and use Julia broadcast as usual:

```julia
first.(coordinates.(centroid.(geotable.geometry)))

```

Is it clearer now how to operate on these geotables?

---

<div class="post-metadata">

**Author:** ![alex-s-gardner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex-s-gardner/32/30210_2.png) [@alex-s-gardner](https://discourse.julialang.org/u/alex-s-gardner)\
**Post date:** [August 25, 2023, 7:13pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/9 "2023-08-25T19:13:26Z")

</div>

Hi @juliohm I don’t see a @transform macro in GeoTables or in Base.

---

<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 25, 2023, 7:29pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/10 "2023-08-25T19:29:25Z")

</div>

That is because GeoTables.jl loads a geospatial data set in the context of the GeoStats.jl framework. If you need the macro for these transformations you need to load the main module of the framework or at least GeoStatsBase.j

---

<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 25, 2023, 7:59pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/11 "2023-08-25T19:59:45Z")

</div>

@chadagreene let me ask you an even better question. What is your ultimate goal trying to get the x-y coordinates? It is often the case that you don’t need to access any coordinates to perform the task you have in mind.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [August 25, 2023, 8:05pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/12 "2023-08-25T20:05:01Z")

</div>

As I mentioned before from mobile, reading vector data with GMT is very handy

```julia
using GMT

D = gmtread("1SDH_201706-Ross1.gpkg")
Attributes: Dict("s1name" => "S1A_EW_GRDM_1SDH_20170622T091354_20170622T091458_017144_01C95D_AB94;S1A_EW_GRDM_1SDH_20170610T091353_20170610T091458_016969_01C402_C103;", "name" => "Ross1", "DATE_" => "201706", "version" => "v1", "updated" => "20221206")
BoundingBox: [-481331.4907714399, -351057.6779158652, -1.2323228635592712e6, -1.1820828635592712e6]
PROJ: +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
WKT: PROJCS["WGS 84 / Antarctic Polar Stereographic",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",-71],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1]]
3288×2 GMTdataset{Float64, 2}
  Row │ X Y
      │ Float64 Float64
──────┼────────────────────────
    1 │ -3.51058e5 -1.18208e6
    2 │ -3.51067e5 -1.18208e6
    3 │ -3.51067e5 -1.18224e6
    4 │ -3.51107e5 -1.18224e6
    5 │ -3.51107e5 -1.18244e6
,,,

```

extracting a subset is straightforward. Get first 5 points

```julia
D[1:5,:]
5×2 Matrix{Float64}:
 -3.51058e5 -1.18208e6
 -3.51067e5 -1.18208e6
 -3.51067e5 -1.18224e6
 -3.51107e5 -1.18224e6
 -3.51107e5 -1.18244e6

```

---

<div class="post-metadata">

**Author:** ![chadagreene](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chadagreene/32/34924_2.png) [@chadagreene](https://discourse.julialang.org/u/chadagreene)\
**Post date:** [August 25, 2023, 8:20pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/13 "2023-08-25T20:20:56Z")

</div>

@juliohm Good question. The answer: I want to do _everything_ with the raw data. For example,

- Whenever I get a new dataset, the first thing I typically want to do is look at the data via `plot(x,y)`.

- Next, I usually want to find out what’s going on at those data points, and that might mean interpolating a gridded field of elevation at every x,y coordinate. In MATLAB, that would look like `z = interp2(x_grid,y_grid,z_grid,x,y)` and then I could do `scatter(x,y,10,z)` to see how elevation changes along the lines.

- I might even want to move the data points around. For example, this dataset represents points on a glacier, and glaciers are constantly moving. So I might want to use an ice velocity field to move the x,y points upstream or downstream.

For most tasks that I want to accomplish, I think it makes the most sense to work directly with the x,y coordinate data.

---

<div class="post-metadata">

**Author:** ![chadagreene](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chadagreene/32/34924_2.png) [@chadagreene](https://discourse.julialang.org/u/chadagreene)\
**Post date:** [August 25, 2023, 8:21pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/14 "2023-08-25T20:21:30Z")

</div>

Thanks @joa-quim, that’s very helpful!

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [August 25, 2023, 8:27pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/15 "2023-08-25T20:27:50Z")

</div>

> [@chadagreene](#):
>
> In MATLAB, that would look like `z = interp2(x_grid,y_grid,z_grid,x,y)`

```julia
using GMT

D2 = grdtrack(D, grid="the_grid_to_innerpolate_from")

```

---

<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 25, 2023, 9:03pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/16 "2023-08-25T21:03:56Z")

</div>

I will answer all your queries bellow without the low-level access to coordinates.

Simply assume that you have loaded the data with:

```julia
gdf = GeoTables.load("fronts/1SDH_201706-Ross1.gpkg")

```

Also assume the following packages loaded:

```julia
using GeoStats # stack for geospatial data science
using GeoTables # package for IO with old GIS formats
import GLMakie as Mke # visualization stack

```

## Q1.

> [@chadagreene](#):
>
> Whenever I get a new dataset, the first thing I typically want to do is look at the data via `plot(x,y)`.

## A1.

```julia
viz(gdf.geometry)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/9/c9a14d7527d90820302718bdfb2721a89ccb2ac0.png)

## Q2.

> [@chadagreene](#):
>
> Next, I usually want to find out what’s going on at those data points, and that might mean interpolating a gridded field of elevation at every x,y coordinate

## A2.

Check the documentation of estimation problems:

[https://juliaearth.github.io/GeoStats.jl/stable/problems.html#Estimation](https://juliaearth.github.io/GeoStats.jl/stable/problems.html#Estimation)

No need to access coordinates explicitly either.

## Q3.

> [@chadagreene](#):
>
> I might even want to move the data points around.

## A3.

```julia
tdf = gdf |> Translate(1000, 1000)

viz(gdf.geometry)
viz!(tdf.geometry, color = :red)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/3/b318286a7c713eaf7997483932d9d4501628da4f.png)

As you can see, you don’t need access to the underlying coordinates if you learn the high-level framework. It is our fault that the framework is not well-understood yet, and we are working to educate more people about it.

In the following months you will see some pretty big improvements on that regard. Stay tuned.

---

<div class="post-metadata">

**Author:** ![alex-s-gardner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex-s-gardner/32/30210_2.png) [@alex-s-gardner](https://discourse.julialang.org/u/alex-s-gardner)\
**Post date:** [August 25, 2023, 9:08pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/17 "2023-08-25T21:08:14Z")

</div>

@juliohm great explanation. One more question: how could one translate individual points?

---

<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 25, 2023, 9:16pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/18 "2023-08-25T21:16:56Z")

</div>

The same Translate transform works with arbitrary collections of geometries.

Let’s take the only geometry in this dataset (a MultiRope) and see what it is made of:

```julia
g = gdf.geometry[1]

v = vertices(g)

```

```julia
3288-element Vector{Point2}:
 Point(-351057.6779158652, -1.1820828635592712e6)
 Point(-351066.50089631113, -1.1820828635592712e6)
 Point(-351066.50089631113, -1.1822428635592712e6)
 Point(-351106.50089631113, -1.1822428635592712e6)
 Point(-351106.50089631113, -1.1824428635592712e6)
 Point(-351146.50089631113, -1.1824428635592712e6)
 Point(-351146.50089631113, -1.1825628635592712e6)
 Point(-351186.50089631113, -1.1825628635592712e6)
 Point(-351186.50089631113, -1.1826028635592712e6)
 Point(-351226.50089631113, -1.1826028635592712e6)
 Point(-351226.50089631113, -1.1826828635592712e6)
 Point(-351266.50089631113, -1.1826828635592712e6)
 Point(-351266.50089631113, -1.1827228635592712e6)
 Point(-351306.50089631113, -1.1827228635592712e6)
 Point(-351306.50089631113, -1.1828028635592712e6)
 Point(-351346.50089631113, -1.1828028635592712e6)
 Point(-351346.50089631113, -1.1828428635592712e6)
 Point(-351386.50089631113, -1.1828428635592712e6)
 ⋮
 Point(-480466.50089631113, -1.2312428635592712e6)
 Point(-480506.50089631113, -1.2312428635592712e6)
 Point(-480506.50089631113, -1.2312028635592712e6)
 Point(-480546.50089631113, -1.2312028635592712e6)
 Point(-480546.50089631113, -1.2311628635592712e6)
 Point(-480586.50089631113, -1.2311628635592712e6)
 Point(-480586.50089631113, -1.2310828635592712e6)
 Point(-480626.50089631113, -1.2310828635592712e6)
 Point(-480626.50089631113, -1.2310428635592712e6)
 Point(-480666.50089631113, -1.2310428635592712e6)
 Point(-480666.50089631113, -1.2309628635592712e6)
 Point(-480706.50089631113, -1.2309628635592712e6)
 Point(-480706.50089631113, -1.2309228635592712e6)
 Point(-480786.50089631113, -1.2309228635592712e6)
 Point(-480786.50089631113, -1.2308828635592712e6)
 Point(-481306.50089631113, -1.2308828635592712e6)
 Point(-481306.50089631113, -1.2309228635592712e6)
 Point(-481331.4907714399, -1.2309228635592712e6)

```

As you can see the single geometry of this dataset (a MultiRope) is made of 3k points. You can transform geometries, points by placing them into a geospatial domain for the transformation.

In this case we can simply place all points into a `PointSet`:

```julia
newpointset = PointSet(v) |> Translate(1000, 1000)

```

Notice that translate is one of more than 25 transforms defined in the framework:

[https://juliaearth.github.io/GeoStats.jl/stable/transforms.html](https://juliaearth.github.io/GeoStats.jl/stable/transforms.html)

We could also add methods that take raw vectors of geometries in these transforms, we just didn’t add because most of the time we are operating on entire domains with smart optimizations.

---

<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:** [February 28, 2024, 12:50pm UTC](https://discourse.julialang.org/t/accessing-x-y-values-in-a-geopackage/103196/19 "2024-02-28T12:50:59Z")

</div>

Update for those coming late here:

Please replace `GeoTables.load` by `GeoIO.load`. The function has been migrated to the GeoIO.jl package a long time ago.