Julia programming: My Journey Begins at EGU

The Beginning

It's been a long time since I first heard about Julia. Despite my interest, I never got around to actually studying it. That changed recently when I attended a short course on Julia at the last EGU General Assembly. Inspired by what I learned, I returned home determined to dive into programming with Julia.

First Steps

I started with a basic book on Julia, aiming to build a solid foundation. The simplicity and power of the language fascinated me, and soon, I felt ready to take a more hands-on approach.

My Project: A Julia Package

As I gained confidence, I began writing my own package. This project is still a work in progress on github. The package includes:

• Reading Time Series: It can read a time series from a CSV file, handling special formatting and metadata.
• Goodness of Fit Functions.
• Graph Operations on river networks.(in progres)
• Geometry/Geography Operations.

Focus on Geography

In this blog post, I want to delve into the last point — geometry and geography operations. My initial goal was simple: create a grid of points within a specified region. This function lays the groundwork for more complex spatial analyses.

Essential Libraries for Geospatial Processing in Julia

• ArchGDAL: This library is crucial for reading and manipulating raster and vector data, forming the backbone of our data handling.
• GeoDataFrames: Built upon ArchGDAL, GeoDataFrames is perfect for reading shapefiles, allowing for easy integration and manipulation of vector data. It also allows the reprojection of features. In this case, pay attention to use the order argument properly, as you can see in this thread.
• GeoFormatTypes: This library manages specific geospatial data types, such as coordinate reference systems (CRS). It's a utility library.
• LibGEOS: For performing complex spatial operations, LibGEOS provides a robust set of tools that are indispensable for any geospatial analysis, particularly when constructing and manipulating geometries, such as unions, intersections, overlaps, etc.

Here is an example of how to reproject a DataFrame with geometry columns (obtained by reading a shapefile with GeoDataFrames) using ArchGDAL and GeoDataFrames, plus GeoFormatTypes:

using GEOFRAMEjl
import ArchGDAL as GDAL
import GeoDataFrames as GDF
import GeoFormatTypes as GFT
using DataFrames

# read the vector file and create a DataFrame with a geometry column
basins_file = ArchGDAL.read("./basin_km3.shp")
basins_layer = ArchGDAL.getlayer(basins_file)
source_crs = ArchGDAL.getspatialref(basin_layer)
target_crs = ArchGDAL.importEPSG(32632, order=:trad)

ArchGDAL.createcoordtrans(source_crs, target_crs) do transform
    for i in 1:nrow(basins)
        geom = basins[i, 1]
        ArchGDAL.transform!(geom, transform)
    end
end

# with GeoDataFrames
basins = GDF.read("./basin_km3.shp")

function get_epsg(df)
    a = GDF.metadata(df)["crs"]
    return convert(GFT.EPSG, a).val[1]
end

source_crs = GFT.EPSG(get_epsg(basins))
target_crs = GFT.EPSG(32632)

basins.geometry = GDF.reproject(basins.geometry, source_crs, target_crs, order=:trad)

Creating a Grid with Julia

Now, let's dive into how I utilize these tools to create a grid of points within a specified region:

using GEOFRAMEjl
import GeoDataFrames as GDF
import ArchGDAL as GDAL
using GeoArrays

# create a squared grid
step = 500
full_grid = Geo.create_grid(614000, 680500, step, 5112000, 5159000, step)

# read polygons (3km2)
basins = GDF.read("./basin_km3.shp")
grid = Geo.filter_points(full_grid, basins, buffer=b)
dem = GeoArrays.read("./dtm_tot_5.tif", band=1)
grid = Geo.get_value_from_raster(grid, dem)
Geo.save_grid(grid, "./grid.shp", 32632)

This generates:

For reprojection, pay attention to the order; more details are available in the thread linked above.