Merge (mosaic) of rasters changes resolution

Question

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:

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:

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

I suspect that it's your call to `resample()` Unfortunately, I can't know for sure because your code isn't reproducible as you have posted it. Would you mind posting the output of `print(ResampleReference_Raster)`, and `print(SDSasRaster)`. This will show the resolution and extent of the rasters. — Sean McKenzie, Jul 23 '21 at 15:38
On the subject of projections, if you are trying to reproject a raster, I believe the correct function is `projectRaster()`. As a best practice, you may want to consider defining your projection strings with their EPSG codes using the syntax `crs(SDS)<-crs("+init=epsg:4326")` for WGS84 — Sean McKenzie, Jul 23 '21 at 15:44
I added the prints to the post, ResampleReference_Raster and SDSasRaster have different resolutions, because of that I used resample. The projection seems correct, the only difference is the inclination of the pixels in Panoply. — Just_4n0th3r_Pr0gr4mm3r, Jul 23 '21 at 15:57
I think the problem is the extent. The reference extent that I am using ignores some islands. I need the rasters to have the same extent because I want to stack them and do the mosaic, how can I obtain a extent that include both tiles? — Just_4n0th3r_Pr0gr4mm3r, Jul 23 '21 at 17:17

score 0 · Accepted Answer · answered Aug 04 '21 at 21:37

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)

Merge (mosaic) of rasters changes resolution

1 Answers1