Downscaling sentinel-2 bands to 10m using R language

Question

I'm trying to calculate NDRE using sentinel-2 bands in R language.

The formula for NDRE = (nir-re)/(nir+re)
nir- Near InfraRed (Band8)
re - RedEdge (Band5)

My Code:

library(raster)
library(RStoolbox)
re_path <- "D:/R/T43PHS_20190223T050811_B05.jp2"
nir_band <- "D:/R/T43PHS_20190223T050811_B08.jp2"
re <- raster(re_path)
nir <- raster(nir_band)
plot((nir-re)/(nir+re), main="NDRE")
writeRaster(x = ((nir-re)/(nir+re)),
            filename="D:/R/T43PHS_20190223T050811.tif",
            format = "GTiff", # save as a CDF
            datatype='FLT4S'
)

But there seems to be an error due to difference in Bands5 and Band8 resolution.

Error in compareRaster(e1, e2, extent = FALSE, rowcol = FALSE, crs = TRUE, : different resolution

You can download Band5 and Band8 Here

I want to convert or downscale the 20m band into 10m band using R language and then calculate the indices, I tried using resample() in R I got the output "tiff" file but there is so much loss of information.

Thank you in advance

wouldn't the other way make more sense? getting 10m raster out of a 20m resolution seems to be a questionable endeavor. why not scale the 10m resolution to a 20m resolution? or is that too imprecise? — D.J, Mar 15 '21 at 07:48
of course we can do that but I need accurate data, for accurate data downsizing 20m to 10m seems to be precise. — Repsol Rider, Mar 15 '21 at 10:37
Now there is the `sen2r` package that can be used to download, convert, resample, and calculate indexes of Sentinel-2 image. it is done quite well and once you have mastered the workflow it becomes easy and intuitive. However you can use band 8 instead of band 8A. Those bands are very similar but the 8 band has a native resolution of 10 m — Elia, Apr 30 '21 at 09:50

Mark Thompson · Answer 1 · 2022-04-28T22:00:28.593

You need to tranlate jp2 into the correct format.

Band Resolutions:

B2:B4,B8 = 10 m
B5:B7 = 20 m

Adapted the "jp2_to_GTiff" function from here, with revision:

jp2_to_GTiff <- function(jp2_path) {
      sen2_GDAL <- raster(readGDAL(jp2_path))
      names(sen2_GDAL) <- as.character(regmatches(jp2_path,gregexpr("B0\\d", jp2_path)))
      return(sen2_GDAL)
    }

## Your files:

re_path <- "D:/R/T43PHS_20190223T050811_B05.jp2"
nir_band <- "D:/R/T43PHS_20190223T050811_B08.jp2"

## Call the function
re <- jp2_to_GTiff(re_path)
nir <- jp2_to_GTiff(nir_band)

## Resample from raster package, then write:
re_10mr <- resample(re, nir, method='bilinear')

writeRaster(x = ((nir-re_10mr)/(nir+re_10mr)),
            filename="D:/R/T43PHS_20190223T050811.tif",
            format = "GTiff", # save as a CDF
            datatype='FLT4S')

Note that the function uses a default GDAL conversion. The command line utility can be used as described here.

score 0 · Answer 2 · answered Sep 03 '22 at 19:19

You can use area-to-point regression Kriging (ATPRK) to downscale S2 using the atakrig package. ATPRK consists of two steps:

regression
area-to-point Kriging on the regression's residuals

Assuming your coarse and fine bands are the red and nir, respectively, and the CRS is in meters, then:

#linear regression
library(raster)

#extract prediction part of a lm model as a raster
red = raster ("path/red.tif")
vals_red <- as.data.frame(values(red))
red_coords = as.data.frame(xyFromCell(red, 1:ncell(red))
combine <- as.data.frame(cbind(red_coords, vals_red))

nir = raster ("path/nir.tif")
nir <- resample(nir, red, method="bilinear")
vals_nir <- as.data.frame(values(nir))

s = stack(red, nir)

block.data <- as.data.frame(cbind(combine, vals_nir))
names(block.data)[3] <- "red"
names(block.data)[4] <- "nir"

block.data <- na.omit(block.data)

block.data = read.csv("path/block.data.csv")

model <- lm(formula = red ~ nir, data = block.data)
#predict to a raster
r1 <- predict(s, model, progress = 'text')
plot(r1)
writeRaster(r1, filename = "path/lm_pred.tif")

#extract residuals
map.resids <- as.data.frame(model$residuals)
x = as.data.frame(block.data$x)
y = as.data.frame(block.data$y)
map.resids <- SpatialPointsDataFrame(data=map.resids, coords=cbind(x,y)) 
gridded(map.resids) <- TRUE
r <- raster(map.resids)
raster::crs(r) <- "EPSG:...." #your EPSG in meters
writeRaster(r, 
            filename = "path/lm_resids.tif", 
            format = "GTiff", 
            overwrite = T)

#area-to-point Kriging
library(atakrig)

block.data = read.csv("path/coords.csv")#csv with 2 columns (x,y) taken from the high resolution raster

nir = raster("path/nir.tif")
rsds = raster("path/lm_resids.tif")

#raster to points
nir.d = discretizeRaster(nir , 10, type = "value")
rsds.d = discretizeRaster(rsds, 10, type = "value")

# point-scale cross-variogram
aod.list <-list(rsds = rsds.d, nir = nir .d)
sv.ck <-deconvPointVgmForCoKriging(aod.list, 
                                   model = "Sph",
                                   rd = 0.5,
                                   maxIter = 50,
                                   nopar = F) #play with the rd and model

ataStartCluster()
pred.atpok <- atpKriging(rsds.d, 
                         block.data, 
                         sv.ck$nir.rsds, 
                         showProgress = T,
                         nopar = F)
ataStopCluster()

# convert result to raster for atp
pred.atpok.r <-rasterFromXYZ(pred.atpok[,-1])
#plot(pred.ataok.r)
raster::crs(pred.atpok.r) <- "EPSG:...." #your EPSG (in meters)
writeRaster(pred.atpok.r$pred, 
            'path/atpk.tif', overwrite = T)


pred = raster('path/lm_pred.tif')
atpk = raster('path/atpk.tif')
pred = resample(pred, atpk, method = 'bilinear')

red_atprk = pred + atprk

red_atprk[red_atprk <= 0] <- 0

writeRaster(red_atprk,
            filename = "path/red_atprk.tif")

Important note: In order ATPRK to produce good results, your bands (coarse and fine) needs to be highly correlated. You can check that by investigating the results of the lm (i.e., high R2 and low AIC or BIC).

Downscaling sentinel-2 bands to 10m using R language

2 Answers2