# GeoArrays band math

**URL:** <https://discourse.julialang.org/t/geoarrays-band-math/65947>\
**Category:** Geo\
**Created:** [August 6, 2021, 3:43pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947 "2021-08-06T15:43:44Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![KrOstir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/krostir/32/27919_2.png) [@KrOstir](https://discourse.julialang.org/u/KrOstir)\
**Post date:** [August 6, 2021, 3:43pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/1 "2021-08-06T15:43:44Z")

</div>

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:

```julia
# 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?

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [August 6, 2021, 4:18pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/2 "2021-08-06T16:18:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![KrOstir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/krostir/32/27919_2.png) [@KrOstir](https://discourse.julialang.org/u/KrOstir)\
**Post date:** [August 6, 2021, 4:35pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/3 "2021-08-06T16:35:13Z")

</div>

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](https://apps.sentinel-hub.com/eo-browser/)

 ![msedge_210806_182739](https://global.discourse-cdn.com/julialang/original/3X/e/3/e3613918236f9e30f437e25b1e02e6339aac726f.jpeg)

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

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [August 6, 2021, 6:01pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/4 "2021-08-06T18:01:32Z")

</div>

Didn’t know this one. Thanks.

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [August 6, 2021, 6:31pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/5 "2021-08-06T18:31:24Z")

</div>

> [@KrOstir](#):
>
> This works, but I have a feeling that this could be easier (and more Julian).

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

```julia
# 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.

> [@KrOstir](#):
>
> Can I extract bands (slice) from the array as an GeoArray?

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

```julia
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.

> [@KrOstir](#):
>
> Is it possible to do the band maths directly on the GeoArray?

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

cc @evetion

---

<div class="post-metadata">

**Author:** ![KrOstir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/krostir/32/27919_2.png) [@KrOstir](https://discourse.julialang.org/u/KrOstir)\
**Post date:** [August 8, 2021, 8:55am UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/6 "2021-08-08T08:55:33Z")

</div>

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

```julia
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)`?

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [August 8, 2021, 9:21am UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/7 "2021-08-08T09:21:04Z")

</div>

> [@KrOstir](#):
>
> As I understand, `ArchGDAL.getband` should be extended to `GeoArray.readband(fn, band)` ?

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.

---

<div class="post-metadata">

**Author:** ![KrOstir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/krostir/32/27919_2.png) [@KrOstir](https://discourse.julialang.org/u/KrOstir)\
**Post date:** [August 13, 2021, 3:15pm UTC](https://discourse.julialang.org/t/geoarrays-band-math/65947/8 "2021-08-13T15:15:20Z")

</div>

Issue raised.  
[Read single band from raster · Issue #72 · evetion/GeoArrays.jl (github.com)](https://github.com/evetion/GeoArrays.jl/issues/72)