Julia for Geospatial part 3 - reading, displaying and writing raster data

This is part three of the Julia for Geospatial series. In part one, available here, I introduced Julia and showed you how to install it, plus how to get the packages you will need to work with raster data. In part two, available here, I directly compared the speed of Python and Julia by opening a raster data and reading it directly to an array. Now we are going to look at reading, displaying and writing raster data using Julia.

There is an associated Jupyter Notebook associated to this blog. Please go ahead and download it. I have added more code and comments to it than in this post.

As we go through the following sections, I'll introduce new concepts and, in part, refactor the code. If there are any questions please let me know. This is a new topic to me as well, but I have found my Python knowledge of GDAL to be very useful. I hope you do too.

Reading an image

We looked at this partly last time. Let’s quickly refresh. The code below will open a raster called Isle_wight.tif. Select the first band and read it to an array called new1.

using ArchGDAL; const AG = ArchGDAL

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

We can get other useful properties from the raster, which we will use later, like getting the projection and geotransformation parameters:

ref = AG.getproj(dataset)
geotransform = AG.getgeotransform(dataset)

Getting the width and height of the dataset:

width = ArchGDAL.width(dataset)
height = ArchGDAL.height(dataset)

And getting the number of raster bands in a dataset (remembering that we often can have several bands of data).

number_rasters = (ArchGDAL.nraster(dataset))

Hopefully, if you are familiar with GDAL, then this will make sense. There are more useful bits in this tutorial: https://github.com/JuliaGeo/GDAL.jl/blob/master/test/tutorial_raster.jl

Plotting Rasters

To plot rasters we need the equivalent of matplotlib.pyplot.imshow(). I found that the colors Package (installed in part one) allows us to do this. To plot a three band RGB image in Julia, use the following code:

using Colors
ArchGDAL.registerdrivers() do    
ArchGDAL.read("isle_wight.tif") do dataset
    number_rasters = (ArchGDAL.nraster(dataset))
    width = ArchGDAL.width(dataset)
    height = ArchGDAL.height(dataset)

    m = Array{Float64}(undef, height, width, number_rasters)

    
    for i = 1:number_rasters
      band = ArchGDAL.getband(dataset,i)
        b = ArchGDAL.read(band)/255
        b_n = permutedims(b, (2,1))
        m[:, :, i] = b_n
    end
    
    imgc2 = colorview(RGB, permuteddimsview(m, (3,1,2)))

    end
end

There is a fair bit going on here. We are reading in the Isle_wight.tif as a dataset and extracting its width and height. We are then using this information to build an array (m) of these dimensions as a float64 type.

Now that we have an array filled with undefined values, we use a for loop to iterate over the bands in the image, and convert each one to a float (b) - this is a plotting issue with the Colors package, it seems. We then change the dimensions of that array (b_n) so it displays correctly (see the accompanying notebook for a more detailed explanation) and then finally write this band to the array m.

Finally, we plot (imgc2) using the colorview method, assigning our data as a RGB. All being well this should be your output:

Writing raster data

Now that we have read and displayed this RGB data, let’s use ArchGDAL to write out a single band to a new tif file. The code is below:

using ArchGDAL; const AG = ArchGDAL ## make AG the shortcut to ArchGDAL, by doing this I call AG. rather than ArchGDAL.
AG.registerdrivers() do
    AG.read("Isle_wight.tif") do dataset
    band1 = ArchGDAL.getband(dataset,1) ## read in 1 band
    ref = AG.getproj(dataset) ## get the project
    new1 = ArchGDAL.read(band1) ## get the data into an array
    geotransform = AG.getgeotransform(dataset) ## get the geotransformation propoerties
        
    raster = AG.unsafe_create(
        "band1.tif",
        AG.getdriver("GTiff"),
        width = ArchGDAL.width(dataset),
        height = ArchGDAL.height(dataset),
        nbands = 1,
        dtype = UInt8
    )
    ## assign the projection and transformation parameters
    AG.setgeotransform!(raster, geotransform)
    AG.setproj!(raster, ref)
    
    ## write the raster    
    AG.write!(
        raster,
        new1,  # image to "burn" into the raster
        1,      # update band 1
    )
    AG.destroy(raster)
    end
end

The first part of this script is similar to before; we read in the data to an array and extract the properties of the raster that we need to write out. Using AG.unsafe_create we parse in the filename to be written along with its parameters (width, height, number of band and type). This is assigned to a variable called raster and it is this variable that we then assign the projection and geotransformation parameters to. This is exactly like we would in Python.

Finally we call AG.write! and parse it our raster (ArchGDAL object), the array containing the data (this could be the result of a calculation, for example, but in our case it’s a copy of band1 of the Isle_wight.tif) and the band. By calling AG.destroy we remove the raster object from memory and end the script.

You could also read and plot this back in to check that the data has been written. It’s probably easiest to use a GIS to do that, and best to use QGIS 3 (of which I have written a couple of books).

This should give you the foundation to go on and do much more with raster datasets.

Recap

In part three we have:

• Read
• Plotted
• Written

raster data using ArchGDAL.

All the code for this series is on my GitHub pages. Please free feel to go and experiment with Julia. I'd love to know how you got on.

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