用R语言将哨兵-2波段降到10米
用R语言将哨兵-2波段降到10米
提问于 2021-03-15 06:49:34
NDRE = (nir-re)/(nir+re)的公式
nir-近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'
)但是，由于difference in Bands5和Band8 分辨率，似乎存在一个错误。
compareRaster中的错误(e1，e2，区段= FALSE，rowcol = FALSE，crs = TRUE，不同分辨率
您可以下载Band5和Band8 这里
我想用R语言将20m波段转换或缩小到10m波段，然后计算指标，在R中尝试使用resample()，得到输出的"tiff“文件，但信息损失很大。
提前谢谢你
回答 2
Stack Overflow用户
发布于 2022-04-28 02:35:51
您需要将jp2转换成正确的格式。
波段分辨率：
- B2:B4，B8 = 10米
- B5:B7 = 20米
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转换。命令行实用程序可以像描述的这里那样使用。
Stack Overflow用户
发布于 2022-09-03 19:19:18
您可以使用面积对点回归克里格(ATPRK)缩小规模的S2使用atakrig软件包.该方案由两个步骤组成：
- 回归
- 回归残差的面积对点Kriging
假设您的粗带和细带分别是红色和nir，并且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")重要注意事项：为了取得好的效果，你的波段(粗的和细的)需要高度相关的。您可以通过调查lm的结果(即高R2和低AIC或BIC)来检查这一点。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/66633464
复制相关文章
相似问题
腾讯云开发者
