Уменьшение масштаба диапазона sentinel-2 до 10 м с использованием языка R

Question

Я пытаюсь рассчитать NDRE, используя группы sentinel-2 на языке R.

Формула для NDRE = (nir-re) / (nir+re)
nir - Near InfraRed (Band8)
re - RedEdge (Band5)

Мой код:

      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'
)

Но , кажется, ошибка из - за разницы в Bands5 и Band8 разрешении .

Ошибка в compareRaster(e1, e2, extract = FALSE, rowcol = FALSE, crs = TRUE,: другое разрешение

Вы можете скачать Band5 и Band8 здесь

Я хочу преобразовать или уменьшить диапазон 20 м до диапазона 10 м с помощью языка R, а затем вычислить индексы, я попытался использовать resample() в R я получил выходной файл "tiff", но там очень много информации.

заранее спасибо

Answer 1

Вы можете использовать регрессионный кригинг «площадь-точка» (ATPRK) для уменьшения масштаба S2 с помощьюatakrigупаковка. АТПРК состоит из двух ступеней:

1. регрессия
2. межточечный кригинг по остаткам регрессии

Предполагая, что ваши грубые и тонкие полосы - это красный и нир соответственно, а CRS - в метрах, тогда:

      #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")

Важное примечание : чтобы ATPRK давал хорошие результаты, ваши полосы (грубая и тонкая) должны быть сильно коррелированы . Вы можете убедиться в этом, изучив результатыlm(т.е. высокаяR2и низкийAICилиBIC).

Answer 2

Вам нужно перевести jp2 в правильный формат.

Разрешение полосы:

• В2:В4,В8 = 10 м
• В5:В7 = 20 м

Адаптировал функцию " jp2_to_GTiff " отсюда , с доработкой:

      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')

Обратите внимание, что функция использует преобразование GDAL по умолчанию. Утилиту командной строки можно использовать, как описано здесь .