Search code examples
rmergerasterresolutionmosaic

Merge (mosaic) of rasters changes resolution

I'm merging two MODIS DSR tiles using a R script that I developed, these are the products: https://drive.google.com/drive/folders/1RG3JkXlbaotBax-h5lEMT7lEn-ObwWsD?usp=sharing

So, I open both products (tile h15v05 and tile h16v05) from same date (2019180), then I open each SDS and merge them together (00h from h15v05 with 00h from h16v05 and so on...)

Visualisation on Panoply (using the merge option) of the two products:

enter image description here

Purple square is the location of the division line that separates the two tiles.

With my code I obtain a plot with pixels with different resolution (and different min/max values) and I don't understand why:

enter image description here

I suspect that the results obtained are due to:

1- Changing from Sinusoidal CRS to longlat WGS84 CRS;

2- Using resample (method ngb) to work with mosaic.

My code is extensive, but here are some parts of it:

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

Does anyone have an idea of what could be the cause of this mismatch between results?

print(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)

print(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)

print(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)

Also, some of the islands were ignored by the script (bottom right)...

enter image description here

Solution

  • sorry I didn't reply earlier. So I think you're right that this issue is extent to which you are resampling. I think you might be able to get around this by creating a dummy raster that has the extent of the raster you want to resample, but has the resolution of the raster you want to mosaic 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)