CODEWITHSUNDEEP
rrastersatellite-imagesentinel2sentinelsat
Repsol Rider
Repsol Rider

Reputation: 21

Downscaling sentinel-2 bands to 10m using R language

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

Upvotes: 2

Views: 1196

Answers (2)

Nikos
Nikos

Reputation: 482

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

  1. regression
  2. 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).

Upvotes: 0

Mark Thompson
Mark Thompson

Reputation: 205

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.

Upvotes: 0

Related Questions