栅格的合并(马赛克)改变分辨率
Merge (mosaic) of rasters changes resolution
我正在使用我开发的 R 脚本合并两个 MODIS DSR 瓦片,这些是产品:
https://drive.google.com/drive/folders/1RG3JkXlbaotBax-h5lEMT7lEn-ObwWsD?usp=sharing
所以,我从同一日期 (2019180) 打开两个产品(tile h15v05 和 tile h16v05),然后我打开每个 SDS 并将它们合并在一起(h15v05 的 00h 和 h16v05 的 00h 等等...)
两种产品的 Panoply 可视化(使用合并选项):
紫色方块是分隔两个方块的分割线的位置。
通过我的代码,我获得了一个具有不同分辨率(和不同 min/max 值)像素的图,但我不明白为什么:
我怀疑得到的结果是由于:
1- 从正弦 CRS 更改为 longlat WGS84 CRS;
2- 使用重采样(方法 ngb)处理马赛克。
我的代码很广泛,但这里是其中的一些部分:
# Open scientific dataset as raster
SDSs <- sds(HDFfile)
SDS <- SDSs[SDSnumber]
crs(SDS) <- crs("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
SDSreprojected <- project(SDS, DesiredCRS)
SDSasRaster <- as(SDSreprojected, "Raster")
# Resample SDS based on a reference SDS (SDS GMT_1200_DSR of a first product), I need to do this to be able to use mosaic
SDSresampled <- resample(SDSasRaster,ResampleReference_Raster,method='ngb')
# Create mosaic of same SDS, but first convert stack to list to use mosaic
ListWith_SameSDS_OfGroupFiles <- as.list(StackWith_SameSDS_OfGroupFiles)
ListWith_SameSDS_OfGroupFiles.mosaicargs <- ListWith_SameSDS_OfGroupFiles
ListWith_SameSDS_OfGroupFiles.mosaicargs$fun <- mean
SDSmosaic <- do.call(mosaic, ListWith_SameSDS_OfGroupFiles.mosaicargs)
# Save SDSs mosaic stack to netCDF
writeRaster(StackWith_AllMosaicSDSs_OfGroupFiles, NetCDFpath, overwrite=TRUE, format="CDF", varname= "DSR", varunit="w/m2", longname="Downward Shortwave Radiation", xname="Longitude", yname="Latitude", zname="TimeGMT", zunit="GMT")
有谁知道导致结果不匹配的原因是什么?
打印(ResampleReference_Raster)
class : RasterLayer
dimensions : 1441, 897, 1292577 (nrow, ncol, ncell)
resolution : 0.01791556, 0.006942043 (x, y)
extent : -39.16222, -23.09196, 29.99652, 40 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : MCD18A1.A2019180.h15v05.061.2020343034815
values : 227.5543, 970.2346 (min, max)
打印(SDSasRaster)
class : RasterLayer
dimensions : 1399, 961, 1344439 (nrow, ncol, ncell)
resolution : 0.01515284, 0.007149989 (x, y)
extent : -26.10815, -11.54627, 29.99717, 40 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : MCD18A1.A2019180.h16v05.061.2020343040755
values : 0, 0 (min, max)
打印(SDSmosaic)
class : RasterLayer
dimensions : 1441, 897, 1292577 (nrow, ncol, ncell)
resolution : 0.01791556, 0.006942043 (x, y)
extent : -39.16222, -23.09196, 29.99652, 40 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : layer
values : 0, 62.7663 (min, max)
此外,脚本忽略了一些岛屿(右下角)...
抱歉没有早点回复。所以我认为你是对的,这个问题是你重新采样的程度。我认为您可以通过创建一个虚拟栅格来解决这个问题,该虚拟栅格具有您想要重新采样的栅格的范围,但具有您想要镶嵌的栅格的分辨率 to.Try:
dummy<-raster(ext = SDSasRaster@extent, resolution=ResampledReference_Raster@res, crs=SDSasRaster@crs)
SDS2<-resample(SDSasRaster, dummy, method="ngb")
Final<-moasic(SDS2, ResampledReference_Raster, fun=mean)
我正在使用我开发的 R 脚本合并两个 MODIS DSR 瓦片,这些是产品: https://drive.google.com/drive/folders/1RG3JkXlbaotBax-h5lEMT7lEn-ObwWsD?usp=sharing
所以,我从同一日期 (2019180) 打开两个产品(tile h15v05 和 tile h16v05),然后我打开每个 SDS 并将它们合并在一起(h15v05 的 00h 和 h16v05 的 00h 等等...)
两种产品的 Panoply 可视化(使用合并选项):
紫色方块是分隔两个方块的分割线的位置。
通过我的代码,我获得了一个具有不同分辨率(和不同 min/max 值)像素的图,但我不明白为什么:
我怀疑得到的结果是由于:
1- 从正弦 CRS 更改为 longlat WGS84 CRS;
2- 使用重采样(方法 ngb)处理马赛克。
我的代码很广泛,但这里是其中的一些部分:
# Open scientific dataset as raster
SDSs <- sds(HDFfile)
SDS <- SDSs[SDSnumber]
crs(SDS) <- crs("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
SDSreprojected <- project(SDS, DesiredCRS)
SDSasRaster <- as(SDSreprojected, "Raster")
# Resample SDS based on a reference SDS (SDS GMT_1200_DSR of a first product), I need to do this to be able to use mosaic
SDSresampled <- resample(SDSasRaster,ResampleReference_Raster,method='ngb')
# Create mosaic of same SDS, but first convert stack to list to use mosaic
ListWith_SameSDS_OfGroupFiles <- as.list(StackWith_SameSDS_OfGroupFiles)
ListWith_SameSDS_OfGroupFiles.mosaicargs <- ListWith_SameSDS_OfGroupFiles
ListWith_SameSDS_OfGroupFiles.mosaicargs$fun <- mean
SDSmosaic <- do.call(mosaic, ListWith_SameSDS_OfGroupFiles.mosaicargs)
# Save SDSs mosaic stack to netCDF
writeRaster(StackWith_AllMosaicSDSs_OfGroupFiles, NetCDFpath, overwrite=TRUE, format="CDF", varname= "DSR", varunit="w/m2", longname="Downward Shortwave Radiation", xname="Longitude", yname="Latitude", zname="TimeGMT", zunit="GMT")
有谁知道导致结果不匹配的原因是什么?
打印(ResampleReference_Raster)
class : RasterLayer
dimensions : 1441, 897, 1292577 (nrow, ncol, ncell)
resolution : 0.01791556, 0.006942043 (x, y)
extent : -39.16222, -23.09196, 29.99652, 40 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : MCD18A1.A2019180.h15v05.061.2020343034815
values : 227.5543, 970.2346 (min, max)
打印(SDSasRaster)
class : RasterLayer
dimensions : 1399, 961, 1344439 (nrow, ncol, ncell)
resolution : 0.01515284, 0.007149989 (x, y)
extent : -26.10815, -11.54627, 29.99717, 40 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : MCD18A1.A2019180.h16v05.061.2020343040755
values : 0, 0 (min, max)
打印(SDSmosaic)
class : RasterLayer
dimensions : 1441, 897, 1292577 (nrow, ncol, ncell)
resolution : 0.01791556, 0.006942043 (x, y)
extent : -39.16222, -23.09196, 29.99652, 40 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : layer
values : 0, 62.7663 (min, max)
此外,脚本忽略了一些岛屿(右下角)...
抱歉没有早点回复。所以我认为你是对的,这个问题是你重新采样的程度。我认为您可以通过创建一个虚拟栅格来解决这个问题,该虚拟栅格具有您想要重新采样的栅格的范围,但具有您想要镶嵌的栅格的分辨率 to.Try:
dummy<-raster(ext = SDSasRaster@extent, resolution=ResampledReference_Raster@res, crs=SDSasRaster@crs)
SDS2<-resample(SDSasRaster, dummy, method="ngb")
Final<-moasic(SDS2, ResampledReference_Raster, fun=mean)