GeoArrays band math

Specific Domains Geo
1

I am working with an GeoArray and would like to perform band maths with it, i.e. compute NDVI from bands. I am currently doing the following:

# Read Sentinel-2 data
fn = download("https://bitbucket.org/earthobservation/julia-remote-sensing/raw/ccea113bf43ba71d4d1c6f02e060fd2923c982b2/data/KR_S2B_MSIL2A_20210614_20.tif")
geoarray = GeoArrays.read(fn)

# Compute NDVI
b_nir = Float64.(geoarray.A[:,:,7]);
b_r = Float64.(geoarray.A[:,:,3]);
ndvi = (b_nir - b_r) ./ (b_nir + b_r);
ndvi_ga = GeoArray(ndvi)
ndvi_ga.f = geoarray.f
ndvi_ga.crs = geoarray.crs

This works, but I have a feeling that this could be easier (and more Julian). Is it possible to do the band maths directly on the GeoArray? Can I extract bands (slice) from the array as an GeoArray?

2

Sorry, not answering to your question, but how do we get those small GeoTIFF Sentinel-2 data sets? All I find in the Copernicus site are those big (large areas) .jp2 with sizes > 10000 x 10000 pixels.

3

The easiest way is to use EO Browser, zoom to area and download analytical data as GeoTIFF. You have to be registered and logged in to use the download feature.

EO Browser

msedge_210806_182739
msedge_210806_1827391920×972 140 KB

The image in the link has been downloaded from Copernicus, cropped and processed in SNAP. :slight_smile:

2 Likes
4

Didn’t know this one. Thanks.

5

Without making any changes to GeoArray, I can only make it slightly simpler:

# Compute NDVI
b_nir = geoarray.A[:,:,7];
b_r = geoarray.A[:,:,3];
ndvi = (b_nir - b_r) ./ (b_nir + b_r);
ndvi_ga = GeoArray(ndvi, geoarray.f, geoarray.crs)

The division converts it to a Float64 already, and we can use the constructor to set the f and crs fields directly.

This drops the GeoArray type, and makes a copy of the data:

julia> geoarray[:,:,7]
900×900 Matrix{UInt16}

It makes some sense that the AffineMap (f) is lost here, since if a subset of rows or columns are taken, it would be invalid. Probably it’s worth raising an issue about adding a method to ArchGDAL.getband, that would preserve the GeoArray, perhaps as a view.

It is, since these operations are implemented: https://github.com/evetion/GeoArrays.jl/blob/master/src/operations.jl. This probably can be expanded though.

cc @evetion

2 Likes
6

With your help, I have modified the code to work nicely on GeoArrays.

ga_b_nir = GeoArray(Float64.(geoarray.A[:,:,7]), geoarray.f, geoarray.crs);
ga_b_r = GeoArray(Float64.(geoarray.A[:,:,3]), geoarray.f, geoarray.crs);
ga_ndvi = (ga_b_nir - ga_b_r) / (ga_b_nir + ga_b_r);

As you said, division converts to Float64. Since both bands are UInt16, the subtraction (can be negative) causes problems, and thus I have made both bands Float64.

Reading only a band to GeoArray would be an excellent option. As I understand, ArchGDAL.getband should be extended to GeoArray.readband(fn, band)?

2 Likes
7

I was thinking perhaps it may be nice if GeoArrays would export ArchGDAL’s getband, and add a method getband(::GeoArray, i::Int)::GeoArray to it. If you’d like something like that, it’s best to raise an issue on the GeoArrays repository.

1 Like
8

Issue raised.
Read single band from raster · Issue #72 · evetion/GeoArrays.jl (github.com)

1 Like