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