CODEWITHSUNDEEP
rgisterrasentinel2
Meirzev
Meirzev

Reputation: 21

Mosaic of multiband Sentinel-2 rasters with Terra R package results in all bands having the same values

I have an issue that is driving me crazy that I just can't figure out for the life of me. I'm processing large Sentinel-2 tiles and attempting to use terra::mosaic() to put them together. I start with two SpatRasters, notice how the min-max values differ for each band:

> sent1
class       : SpatRaster 
dimensions  : 10980, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6400020  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : T09VVD_2022_mean.tif 
names       :   B02,  B03,   B04,        B05,        B06,       B07, ... 
min values  :     0,    0,    29,   379.8125,   523.3125,   385.875, ... 
max values  : 10456, 9464, 15768, 13930.6875, 15912.1250, 16714.938, ... 

> sent2
class       : SpatRaster 
dimensions  : 10980, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6390240, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : T09VVE_2022_mean.tif 
names       :  B02,  B03,   B04,       B05,        B06,      B07, ... 
min values  :    0,  211,     1,   379.125,   331.3125,   342.75, ... 
max values  : 9104, 8448, 10480, 12210.500, 16052.9375, 16814.38, ...

Now see how all bands have the same values:

sent3 <- mosaic(sent1, sent2)

> show(sent3)
class       : SpatRaster 
dimensions  : 20982, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : spat_fubTXtzTOH0dho5_70406.tif 
names       :        B02,        B03,        B04,        B05,        B06,        B07, ... 
min values  :   355.6161,   355.6161,   355.6161,   355.6161,   355.6161,   355.6161, ... 
max values  : 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, ...

When I try to plot it in R or QGIS all the bands look the same which is obviously not the desired output.

I have also tried to add the files to a SpatRasterCollection and mosaic that but it had the same results:

sent.sprc <- sprc(sent1, sent2)

sent4 <- mosaic(sent.sprc)

> show(sent4)
class       : SpatRaster 
dimensions  : 20982, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : spat_hDQlRCrwBUOApYB_49105.tif 
names       :        B02,        B03,        B04,        B05,        B06,        B07, ... 
min values  :   355.6161,   355.6161,   355.6161,   355.6161,   355.6161,   355.6161, ... 
max values  : 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, 12130.0361, ...

However, merge() works and produces the desired output:

sent5 <- merge(sent1, sent2)

> show(sent5)
class       : SpatRaster 
dimensions  : 20982, 10980, 7  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : 399960, 509760, 6290220, 6500040  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
source      : spat_fqpzxaTDk0yKNGb_50974.tif 
names       :   B02,  B03,   B04,       B05,        B06,      B07, ... 
min values  :     0,    0,     1,   379.125,   331.3125,   342.75, ... 
max values  : 10456, 9464, 15768, 13930.688, 16052.9375, 16814.38, ...

Overall I am working with several files that are all quite large so I'm running on an HPC. I transferred two of the files to my local machine to test it out and it's the same problem. I have the current development version of Terra installed and only have 'terra' and 'here' packages loaded. I tried to make a reproducible example but couldn't figure out how to generate dummy data for the multiband raster. If you have any suggestions please let me know! I can supply two of the sentinel-files if that helps but they are >3gb each.

Upvotes: 1

Views: 465

Answers (2)

Meirzev
Meirzev

Reputation: 21

For anyone following this it looks like this was fixed in a package update, haven't tested it myself yet though...

update:

This is now fixed in the current version (1.7.56). I copied the example data from the above github issue.

> library(terra)
terra 1.7.56

> x1 <- rast(xmin=-110, xmax=-60, ymin=40, ymax=70, res=1, vals=1)
> x2 <- rast(xmin=-110, xmax=-60, ymin=40, ymax=70, res=1, vals=2)
> x  <- c(x1, x2)
> 
> y1 <- rast(xmin=-95, xmax=-45, ymax=60, ymin=30, res=1, vals=3)
> y2 <- rast(xmin=-95, xmax=-45, ymax=60, ymin=30, res=1, vals=4)
> y  <- c(y1, y2)
> 
> m1 <- mosaic(x, y)
> m1
class       : SpatRaster 
dimensions  : 40, 65, 2  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : -110, -45, 30, 70  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
source(s)   : memory
names       : lyr.1, lyr.1 
min values  :     1,     2 
max values  :     3,     4

Now the values are different for each layer as opposed to copying the same data for each layer as it was doing before.

Upvotes: 0

UseR10085
UseR10085

Reputation: 8198

You have not provided the data. So, I cannot test my code. You can try the following code:

library(terra)

setwd("...") #Provide the directory containg the the tif files

files <- list.files(".\\", pattern='tif$',
                    full.names = TRUE) 
x <- sprc(files) 
m <- mosaic(x)

Upvotes: 0

Related Questions