On the research of performance tools for spatial analysis we took a chance on exploring Julia for applied raster algebra together with the ArchGDAL package. This tutorial shows the complete procedure to load rasters in Julia, plot them, extract the band arrays, calculate the NDVI index and export the resulting array as a geospatial raster.

For this tutorial you need Julia installed on you computer. Download Julia from this link: julialang.org

Install the Julia kernel for Jupyter following this link: github.com/JuliaLang/IJulia.jl

Tutorial

Code

# Run this once to install the required packages
#using Pkg
#Pkg.add("ArchGDAL")

# load the required packages
using ArchGDAL

# open the nir and red from a Landsat 8 images
redFile = ArchGDAL.read("../Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B4_clip.tif")
nirFile = ArchGDAL.read("../Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B5_clip.tif")

GDAL Dataset (Driver: GTiff/GeoTIFF)
File(s): 
  ../Landsat8/LC08_L1TP_042035_20180603_20180615_01_T1_B5_clip.tif

Dataset (width x height): 2107 x 1338 (pixels)
Number of raster bands: 1
  [GA_ReadOnly] Band 1 (Gray): 2107 x 1338 (UInt16)

# open the band from the raster
redBand = ArchGDAL.getband(redFile, 1)
nirBand = ArchGDAL.getband(nirFile, 1)

[GA_ReadOnly] Band 1 (Gray): 2107 x 1338 (UInt16)
    blocksize: 2107×1, nodata: nothing, units: 1.0px + 0.0
    overviews:

# read the band arrays
redArray = ArchGDAL.read(redBand)
nirArray = ArchGDAL.read(nirFile) 
redArray[10:20,10:20]*1

11×11 Matrix{Int64}:
 7438  7437  7434  7434  7433  7418  7438  7439  7444  7439  7903
 7434  7437  7458  7450  7435  7436  7434  7431  7451  7451  7594
 7449  7455  7446  7440  7432  7432  7427  7464  7466  7455  7541
 7447  7431  7431  7433  7432  7434  7423  7436  7447  7456  7714
 7427  7444  7434  7416  7435  7445  7419  7421  7436  7420  7780
 7427  7432  7447  7426  7435  7438  7434  7427  7431  7456  7579
 7449  7472  7469  7438  7445  7452  7435  7434  7451  7469  7432
 7441  7443  7446  7461  7466  7443  7443  7434  7420  7443  7469
 7419  7440  7449  7423  7448  7447  7440  7419  7443  7453  7640
 7405  7424  7439  7423  7450  7428  7413  7417  7420  7438  7608
 7429  7431  7433  7424  7435  7440  7414  7420  7422  7444  7429

#explore the array type
typeof(redArray)

Matrix{UInt16} (alias for Array{UInt16, 2})

redArrayFloat = Float64.(redArray)
nirArrayFloat = Float64.(nirArray)
redArrayFloat[10:15,10:15]*1

6×6 Matrix{Float64}:
 7438.0  7437.0  7434.0  7434.0  7433.0  7418.0
 7434.0  7437.0  7458.0  7450.0  7435.0  7436.0
 7449.0  7455.0  7446.0  7440.0  7432.0  7432.0
 7447.0  7431.0  7431.0  7433.0  7432.0  7434.0
 7427.0  7444.0  7434.0  7416.0  7435.0  7445.0
 7427.0  7432.0  7447.0  7426.0  7435.0  7438.0

#define a function for the raster algebra
function ndviCal(red,nir)
    #ndviArray = (nir - red)/(nir + red)
    ndviArray = ((nir .- red)./(red .+ nir))[:,:,1]
    return(ndviArray)
end

ndviCal (generic function with 1 method)

@time ndviBand = ndviCal(redArrayFloat,nirArrayFloat)
ndviBand[10:20,10:20]

0.347859 seconds (628.80 k allocations: 74.390 MiB, 14.36% gc time, 79.66% compilation time)





11×11 Matrix{Float64}:
 0.519494  0.521659  0.523935  0.522175  …  0.525935  0.529058  0.492454
 0.524041  0.522305  0.522978  0.52213      0.524171  0.523974  0.51294
 0.521103  0.515893  0.52125   0.522388     0.527633  0.525099  0.514861
 0.52177   0.521198  0.523883  0.525412     0.528313  0.526859  0.502868
 0.523895  0.522239  0.523706  0.523852     0.524583  0.526771  0.498663
 0.519335  0.520856  0.520276  0.523959  …  0.525691  0.526046  0.51262
 0.520039  0.518169  0.519168  0.520516     0.520944  0.524691  0.523986
 0.522186  0.522073  0.521311  0.520994     0.524313  0.527564  0.522671
 0.523139  0.521066  0.521303  0.522068     0.525697  0.525815  0.5077
 0.527562  0.523675  0.51946   0.522575     0.523886  0.526197  0.511117
 0.525546  0.527561  0.521732  0.525956  …  0.526824  0.526839  0.526287

ArchGDAL.create(
    "../Output/NdviJulia.tif",
    driver = ArchGDAL.getdriver("GTiff"),
    width=ArchGDAL.width(redFile),
    height=ArchGDAL.height(redFile),
    nbands=1,
    dtype=Float64
) do dataset
    ArchGDAL.write!(dataset,ndviBand,1)
    ArchGDAL.setgeotransform!(dataset, ArchGDAL.getgeotransform(redFile))
    ArchGDAL.setproj!(dataset, ArchGDAL.getproj(redFile))
end