# \[ANN\] Announcing GeoDataFrames.jl
**URL:** https://discourse.julialang.org/t/ann-announcing-geodataframes-jl/56735
**Category:** Package Announcements
**Tags:** package, announcement, geo
**Created:** [March 8, 2021, 2:35pm UTC](https://discourse.julialang.org/t/ann-announcing-geodataframes-jl/56735 "2021-03-08T14:35:37Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [March 8, 2021, 2:35pm UTC](https://discourse.julialang.org/t/ann-announcing-geodataframes-jl/56735/1 "2021-03-08T14:35:37Z")
I’d like to announce [GeoDataFrames.jl](https://github.com/evetion/GeoDataFrames.jl), a Table interface for reading and writing common vector formats, such as shapefiles and geopackages. It’s quite small, as the majority comes from the work in [ArchGDAL.jl](https://github.com/yeesian/ArchGDAL.jl/) and [GDAL.jl](https://github.com/JuliaGeo/GDAL.jl) on which this package depends.
The main motivation came from the need to quickly read or write some geo-enabled formats. This was already possible with ArchGDAL, but wasn’t as straightforward as it could be. This was visible here ([0,1]) and on Github ([2,3,4]). The package is still quite young and rough, but works for me (albeit simple) production work.
An example:
Writing:
```julia
using GeoDataFrames; const GDF=GeoDataFrames
using DataFrames
coords = zip(rand(10), rand(10))
df = DataFrame(geom=createpoint.(coords), name="test");
GDF.write("test_points.shp", df)
```
Reading:
```julia
df = GDF.read("test_points.shp")
10×2 DataFrame
Row │ name geom
│ String IGeometr…
─────┼───────────────────────────────────────────
1 │ test Geometry: POINT (0.2360926400353…
2 │ test Geometry: POINT (0.2445619453460…
3 │ test Geometry: POINT (0.4504663468371…
4 │ test Geometry: POINT (0.0886989855586…
5 │ test Geometry: POINT (0.0323344938606…
6 │ test Geometry: POINT (0.1574232985696…
7 │ test Geometry: POINT (0.9677152948776…
8 │ test Geometry: POINT (0.0047328946715…
9 │ test Geometry: POINT (0.7389241862917…
10 │ test Geometry: POINT (0.1207370929831…
```
This package may eventually become obsolete as more of the code is upstreamed into ArchGDAL.
[0] [ArchGDAL.getfield() error](https://discourse.julialang.org/t/archgdal-getfield-error/45770)
[1] [How to create a new shapefile containing a few points - #10 by rio](https://discourse.julialang.org/t/how-to-create-a-new-shapefile-containing-a-few-points/43454/10)
[2] [Construct ArchGDAL. IFeatureLayer via named tuples or DataFrame · Issue #152 · yeesian/ArchGDAL.jl · GitHub](https://github.com/yeesian/ArchGDAL.jl/issues/152)
[3] [Documentation [question] writing a table to a layer · Issue #153 · yeesian/ArchGDAL.jl · GitHub](https://github.com/yeesian/ArchGDAL.jl/issues/153)
[4] [Failed to create a Shapefile following the GDAL Vector API tutorial · Issue #147 · yeesian/ArchGDAL.jl · GitHub](https://github.com/yeesian/ArchGDAL.jl/issues/147)
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [March 8, 2021, 3:58pm UTC](https://discourse.julialang.org/t/ann-announcing-geodataframes-jl/56735/2 "2021-03-08T15:58:57Z")
Thank you @evetion for the great work as usual. Just want to point out that the work that is recently happening around Meshes.jl and GeoStats.jl may potentially lead to a cleaner approach with general Tables.jl tables.
Consider a `table` following the `Tables.jl` API:
```julia
julia> using Tables
julia> table = (a=rand(100), b=rand(100)) # a Tables.jl table
(a = [0.307732211237423, 0.803718835997471, 0.08411713448982439, 0.4830139406337839, 0.89972361406729, 0.4838048624166984, 0.6572940355048909, 0.19799338569165803, 0.7805574255598815, 0.8990875618632652 … 0.07242337288301681, 0.7676977890604322, 0.8262285634679147, 0.03258811191692268, 0.5968713113063921, 0.2685941981073232, 0.4662204960862124, 0.5170359078991946, 0.5012109114845784, 0.11183068325228596], b = [0.9965173138643892, 0.11553726662208041, 0.20159410317223725, 0.7313178905382756, 0.557382816267495, 0.2856813840254393, 0.5016940394671034, 0.7547849171693573, 0.1820419942447178, 0.49498420120152087 … 0.003864021720924482, 0.3285419806520167, 0.26864603835206635, 0.23492724473816962, 0.011024989056165335, 0.8322119339786553, 0.5203882014765839, 0.6062219376724753, 0.8817101294035377, 0.6451100662465474])
```
We can georeference this table on various kinds of `Domain` defined in Meshes.jl, e.g. PointSet, CartesianGrid, GeometrySet, etc. For example, we can georeference on a PointSet:
```julia
julia> using GeoStats.jl # v0.23
julia> 𝒟 = georef(table, PointSet(rand(Point3, 100)))
100 GeoData{3,Float64}
variables
└─a (Float64)
└─b (Float64)
domain: 100 PointSet{3,Float64}
```
This object implements the Meshes.jl `Data` trait, and consequently the Tables.jl API:
```julia
julia> first(Tables.rows(𝒟), 3)
3-element Vector{NamedTuple{(:a, :b, :geometry), Tuple{Float64, Float64, Point3}}}:
(a = 0.307732211237423, b = 0.9965173138643892, geometry = Point(0.9489465238556838, 0.6840987940920451, 0.09563305979228809))
(a = 0.803718835997471, b = 0.11553726662208041, geometry = Point(0.055798761161550026, 0.13459614756211802, 0.06163990143241693))
(a = 0.08411713448982439, b = 0.20159410317223725, geometry = Point(0.39791182546602677, 0.47022283912414164, 0.6568206304457964))
```
where an additional column called `geometry` has been generated on the fly with the points in the PointSet. The magic starts when we want to georeference data on large meshes. For example, we can create a large `CartesianGrid` without allocating any memory:
```julia
julia> grid = CartesianGrid(10, 10)
10×10 CartesianGrid{2,Float64}
minimum: Point(0.0, 0.0)
maximum: Point(10.0, 10.0)
spacing: (1.0, 1.0)
```
and in this case the georerence operation is just a lazy wrapper:
```julia
julia> 𝒟 = georef(table, grid)
100 GeoData{2,Float64}
variables
└─a (Float64)
└─b (Float64)
domain: 10×10 CartesianGrid{2,Float64}
```
However, if we want to treat this data set as table we can:
```julia
julia> first(Tables.rows(𝒟), 3)
3-element Vector{NamedTuple{(:a, :b, :geometry), Tuple{Float64, Float64, Quadrangle{2, Float64, Vector{Point2}}}}}:
(a = 0.307732211237423, b = 0.9965173138643892, geometry = Quadrangle(Point(0.0, 0.0), Point(1.0, 0.0), Point(1.0, 1.0), Point(0.0, 1.0)))
(a = 0.803718835997471, b = 0.11553726662208041, geometry = Quadrangle(Point(1.0, 0.0), Point(2.0, 0.0), Point(2.0, 1.0), Point(1.0, 1.0)))
(a = 0.08411713448982439, b = 0.20159410317223725, geometry = Quadrangle(Point(2.0, 0.0), Point(3.0, 0.0), Point(3.0, 1.0), Point(2.0, 1.0)))
```
Notice how the geometry column now has `Quadrangle` from Meshes.jl built on the fly. This is a much needed feature for large-scale geostatistics on unstructured meshes.
So moving forward, if you plan to interop with the GeoStats.jl stack for statistical work, I suggest to make your `GeoDataFrame` type a concrete implementation of the Meshes.jl `Data` trait. This is all that is needed to take advantage of the full stack of geospatial estimation, simulation, learning methods we are building. Additionally, if you implement this trait, you will get plot recipes for free, which will become Makie.jl recipes at some point.
---
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [March 8, 2021, 6:53pm UTC](https://discourse.julialang.org/t/ann-announcing-geodataframes-jl/56735/3 "2021-03-08T18:53:04Z")
Hi Júlio, thanks for the showcase. I really like the pure Julia geometry approach. Do you use robust predicates in Meshes.jl, like [GitHub - lairez/ExactPredicates.jl: Fast and exact geometrical predicates in the Euclidean plane](https://github.com/lairez/ExactPredicates.jl)?
With regards to interfacing, I agree that traits are the way to go. GeoDataFrames makes direct use of (Arch)GDAL geometries, so it would be good to start there. I think that the long overdue rewrite of [GeoInterface](https://github.com/JuliaGeo/GeoInterface.jl/pull/33) is the correct way to go from ArchGDAL geometries to the types `Point` or `Quadrangle` in Meshes and vica-versa. On top of that we can implement the `Data` interface so the attributes come along as well.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [March 8, 2021, 7:13pm UTC](https://discourse.julialang.org/t/ann-announcing-geodataframes-jl/56735/4 "2021-03-08T19:13:19Z")
We have many robust algorithms implemented in Meshes.jl already but some of the predicates are still on the TODO list.
The plan you wrote above is perfect. I will be working on other fronts meanwhile like the MeshIO.jl package to load more meshes directly in other engineering formats not in GIS.