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?
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?
Would you mind posting a link to an example geopackage?
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.
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 = ???
Try
using GMT
D = gmtread(â€śfile.gpkgâ€ť)
And the coordinates (and the rest) will be plainly visible
@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
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?
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
@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.
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
@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.
using GMT
D2 = grdtrack(D, grid="the_grid_to_innerpolate_from")
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
viz(gdf.geometry)
Check the documentation of estimation problems:
https://juliaearth.github.io/GeoStats.jl/stable/problems.html#Estimation
No need to access coordinates explicitly either.
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.
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.
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.