Scale landsat data values

[Nowosad]

I guess because the 'landsat' dataset has (uncorrected) brightness values, so the NDVI turns out to be beyond the valid range? I suggest adding a box to explain that, but I'm not sure what are the values exactly.

I think the best idea is just to rescale the current landsat data (https://www.usgs.gov/faqs/how-do-i-use-a-scale-factor-landsat-level-2-science-products) and replace it with the new version.

[michaeldorman]

Which product does landsat.tif comes from? Depending on the product, the scaling factor is going to be different. The range in landsat.tif is 5252-31961, so I guess it is from Collection-1, but we need to be sure if we're adding the scaling factors to the text. Alternatively, I can download another Landsat image from the same area.

[Nowosad]

Hi @michaeldorman the data is from https://registry.opendata.aws/usgs-landsat/index.html (see https://github.com/Nowosad/spData/blob/17316316b93833fcd4d9f3f640f713b7c6f191c0/data-raw/data-download.R#L90C3-L90C51).

I am pretty sure that it is Collection 2 -- negative values are probably errors (clouds?). I would suggest not doing the scaling in the text, but rather preparing the file beforehand, and just using it after the scaling in the book. (I am still thinking of doing a similar thing in geocompr).

[michaeldorman]

Hi @Nowosad , Please see the following example of rescaling landsat.tif using collection-2 coefficients:

library(stars)
r = read_stars('landsat.tif')
r = r * 0.0000275 - 0.2
red = r[,,,3,drop=TRUE]
nir = r[,,,4,drop=TRUE]
ndvi = (nir-red)/(nir+red)
plot(ndvi)
plot(ndvi < -1 | ndvi > 1, col=c('grey', 'red'))

Substantial parts of the image turn out to be beyond the -1-1 range, unless I'm missing something. I'm not familiar with this situation with Landsat data, will appreciate your advice on how to describe it in the text, thanks.

[Nowosad]

@michaeldorman -- some values are below 0 after scaling -- this usually represent problematic pixels (e.g., clouds). The calculations work after replacing the values below 0 to 0:

library(terra)
#> terra 1.7.67
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast = rast(multi_raster_file)
multi_rast = (multi_rast * 0.0000275) - 0.2
plot(multi_rast)

multi_rast[multi_rast < 0] = 0
ndvi_fun = function(nir, red){
    (nir - red) / (nir + red)
}
ndvi_rast = lapp(multi_rast[[c(4, 3)]], fun = ndvi_fun)
plot(ndvi_rast)

[michaeldorman]

Got it, thank you very much @Nowosad ! Now fixed in https://github.com/geocompx/geocompy/commit/950ad779be9f6f9df20a23ec61eeca39500002c6

[Nowosad]

Great