In order to disrupt existing practices and really develop something new in Julia, we had to make some hard decisions along the way. One of these decisions relates to how we are willing to interface our framework with existing GIS standards and workflows.
On the one hand, we could have followed the path that was followed by other communities such as Python and R, and focus our energy interfacing with well-tested GIS libraries written in C/C++ (e.g., GDAL, GEOS). This is precisely what the JuliaGeo organization has been doing over the years, and it is an important agenda to bring people from other languages that are used to the OGC standards.
On the other hand, we have young geoscientists and first-time programmers who have never studied GIS before, and who really struggle learning the technology as it is today. The widespread emphasis on machine representation and software engineering has created a gap between the developers and the users of GIS software. A typical gap the Julia programming language helps to close.
We decided to limit our interface with existing GIS technology to input and output (IO) of files. This gives users of the framework the chance to
Import geospatial data stored as simple features
Perform geospatial data science with a rich set of tools
Export results to widely used software (e.g., QGIS, ArcGIS)
It creates an ecosystem where users can become contributors and maintainers of the framework, without any knowledge of a second programming language.
3.1 GeoIO.jl
The GeoIO.jl module can load and save geospatial data on disk in a variety of formats, including the most popular formats in GIS (e.g., .shp, .geojson, .kml, .parquet) thanks to various backend packages spread across various Julia organizations. It is designed for users who just want to get their data ready for geospatial data science.
To load a file from disk, we use GeoIO.load:
usingGeoIOgeotable = GeoIO.load("file.shp")
The function automatically selects the backend based on the file extension, converts the simple features into a geospatial domain, and returns a GeoTable.
To save the GeoTable to disk, possibly in a different format, we use GeoIO.save:
GeoIO.save("file.geojson", geotable)
The module fixes inconsistencies between formats whenever possible. For example, the GeoJSON format writes Date columns as String because the JSON format has no date types. The Shapefile format has its own limitations, etc.
Over time, we expect to improve the ecosystem as a whole by highlighting various issues with available standards and backend implementations.
3.2 File formats
Most GIS file formats do not preserve topological information. This means that neighborhood information is lost as soon as geometries are saved to disk. To illustrate this issue, we consider a geotable over a grid:
Other file formats such as .ply and .msh are widely used in computer graphics to save geospatial data over meshes, and preserve topological information:
beethoven = GeoIO.load("data/beethoven.ply")viz(beethoven.geometry)
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227
3.3 Rationale
Now that we have set expectations for our interface with GIS, let’s address an important question that many readers might have coming from other communities:
Do we gain anything by not adhering to programming interfaces?
The answer is an emphatic YES! It means that we have total freedom to innovate and improve the representation of various geometries and geospatial domains with Julia’s amazing type system. To give a simple example, let’s take a look at the Triangle geometry:
t =Triangle((0, 0), (1, 0), (1, 1))
Triangle
├─ Point(x: 0.0 m, y: 0.0 m)
├─ Point(x: 1.0 m, y: 0.0 m)
└─ Point(x: 1.0 m, y: 1.0 m)
If we treated this geometry as a generic polygon represented by a vector of vertices in memory, like it is done in GeoInterface.jl for example, we wouldn’t be able to dispatch optimized code that is only valid for a triangle:
In Julia, the macro @code_llvm shows the underlying code sent to the LLVM compiler. In this case, the code is the single line ret i8 1, which is the instruction to return the constant integer 1.
Notice how the isconvex function is compiled away to the constant1 (i.e. true) when called on the triangle. The code for a generic polygon is much more complicated and requires runtime checks that are too expensive to afford, especially in 3D.
Another reason to not adhere to a generic interface is that we can store information in the geometry types themselves (e.g., coordinate reference system) that is relevant to design advanced scientific visualization features illustrated in the previous chapter, and to dispatch specialized algorithms from geodesic geometry.
Having cleared that up, we will now proceed to the last foundational chapter of the book, which covers the advanced geometric processing features of the framework.
