Accessing x, y values in a geopackage

chadagreene · August 25, 2023, 4:45pm

Basic question: I can load a geopackage with:

using GeoTables 

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

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

alex-s-gardner · August 25, 2023, 4:52pm

Would you mind posting a link to an example geopackage?

chadagreene · August 25, 2023, 4:53pm

@alex-s-gardner https://download.geoservice.dlr.de/icelines/files/Ross1/Ross1_monthly.zip

juliohm · August 25, 2023, 5:01pm

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

geotable.geometry

The entire DataFrame syntax is supported:

https://juliaearth.github.io/GeoStats.jl/stable/quickstart.html

Worth checking the quickstart above.

chadagreene · August 25, 2023, 5:16pm

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?

using GeoTables 

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

joa-quim · August 25, 2023, 5:32pm

Try

using GMT

D = gmtread(“file.gpkg”)

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

alex-s-gardner · August 25, 2023, 5:50pm

@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:

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

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

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

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

There it is

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

let’s index into table.geometry

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

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

juliohm · August 25, 2023, 6:19pm

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:

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

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

@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:

geotable |> viewer

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

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

Is it clearer now how to operate on these geotables?

alex-s-gardner · August 25, 2023, 7:13pm

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

juliohm · August 25, 2023, 7:29pm

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

juliohm · August 25, 2023, 7:59pm

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

joa-quim · August 25, 2023, 8:05pm

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

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

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

chadagreene · August 25, 2023, 8:20pm

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

chadagreene · August 25, 2023, 8:21pm

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

joa-quim · August 25, 2023, 8:27pm

using GMT

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

juliohm · August 25, 2023, 9:03pm

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

Simply assume that you have loaded the data with:

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

Also assume the following packages loaded:

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

Q1.

A1.

viz(gdf.geometry)

Q2.

A2.

Check the documentation of estimation problems:

https://juliaearth.github.io/GeoStats.jl/stable/problems.html#Estimation

No need to access coordinates explicitly either.

Q3.

A3.

tdf = gdf |> Translate(1000, 1000)

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

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.

alex-s-gardner · August 25, 2023, 9:08pm

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

juliohm · August 25, 2023, 9:16pm

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:

g = gdf.geometry[1]

v = vertices(g)

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:

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

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.

juliohm · February 28, 2024, 12:50pm

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.