Julia for Geospatial part 2 - speed test

Last time, I introduced Julia and walked through some basic instructions on how to setup. If you missed it go here. Compared to my first experience of Python, installing Julia was a breeze. In this series I am exploring how suited Julia is for Geospatial programming. This is not to discourage you away from other languages that have better support, but to show the capabilities of Julia.

This time: speed test

The main selling point for Julia, at least for me, is its speed. I wrote last time that I never really had the problem with the speed of Python that others have had. In this post I wanted to re-run my Fastest image Python post and compare the results with Julia. I'll be comparing GDAL in Python, Rasterio in Python and ArchGDAL in Julia.

For this test I am using a clipped out 99% stretched Sentinel 2a 10m true colour RGB image of the Isle of Wight. We are going to open the image in GDAL and read one band as an array. The image is 25.6mb and is shown in QGIS 3.4 below.

Python 3

The code to use GDAL to read an image into an array called image_gdal is below:

from osgeo import gdal
raster_ds = gdal.Open("Isle_wight.tif", gdal.GA_ReadOnly)
image_gdal = raster_ds.GetRasterBand(1).ReadAsArray()

and for Rasterio reading into an array called arr is below:

import rasterio
with rasterio.open("Isle_wight.tif") as r:
    arr = r.read(1)

Python 3 results

I expect the timings to vary slightly on different computers, but I ran this code using the magic function timeit and the results were always comparable. In terms of speed I couldn't really separate them.

I've found the more I have used Rasterio the more I prefer it to Python GDAL. A personal preference.

Julia 1.1.0

The Julia code for ArchGDAL reading into an array called new1 is below:

using ArchGDAL; const AG = ArchGDAL
using BenchmarkTools

function gdal_time()
AG.registerdrivers() do
    AG.read("Isle_wight.tif") do dataset
    band1 = AG.getband(dataset,1)
    new1 = AG.read(band1)
    end
end
end

The closest I could find to the magic function timeit was BenchmarkTools and this is called by the command:

@btime gdal_time()

So, unlike the above, I've had to build a function to perform my speed test. Perhaps there is a better way?

Julia 1.1.0 results

The results showed that Julia using this example was almost twice as fast.

In the time Python took to open the same image and read to an array with GDAL or Rasterio, Julia would have been able to do it twice. So it looks like (from this test) that Julia is approximately two times faster than Python for this process. That was pretty amazing to me, with no hardware upgrades and no GPU involved; there is the potential to really speed up processing.

Recap

In this post I have revisited the fastest image work. We have looked at:

• Reading in a single Satellite band into an array with Python (using GDAL and Rasterio).
• Reading in a single Satellite band into an array with Julia (using ArchGDAL).

From this evidence it appears the claim that Julia is faster than Python holds up for reading raster data.

Now we know that Julia is fast, in part three I'll step through an example of reading, displaying and writing rasters using Julia.

Take a look at the code here

I am a freelancer able to help you with your projects. I offer consultancy, training and writing. I’d be delighted to hear from you.

I have grouped all my previous blogs (technical stuff / tutorials / opinions / ideas) at https://gis.acgeospatial.co.uk.

Feel free to connect or follow me; I am always keen to talk about Earth Observation.

I am @map_andrew on twitter