rioxarray.merge.merge_arrays - bounds argument

[nmaffe]

I am trying to mosaic 6 .tif tiles using rioxarray.merge.merge_arrays. Each tile is 3601 x 3601. The tiles form a mosaic of 2 x 3. The expected result is 10801 x 7201 since these tiles share the same adjacent borders. I have 6 .tif files that I can successfully mosaic to the expected shape 10801 x 7201. I also have other 6 mask .tif tiles of the exact same shape and same boundaries that I'd like to mosaic into a mosaic_mask. I don't understand why in this latter case the size of mosaic_mask is 10802 x 7202 and not 10801 x 7201 as before.

from rioxarray.merge import merge_arrays

folder1 = 'folder1/'
files = ['ASTGTMV003_N45E007_dem.tif', 'ASTGTMV003_N45E008_dem.tif', 'ASTGTMV003_N45E009_dem.tif',
        'ASTGTMV003_N46E007_dem.tif', 'ASTGTMV003_N46E008_dem.tif', 'ASTGTMV003_N46E009_dem.tif']

src_files_to_mosaic = []
for file in files:
    src = rioxarray.open_rasterio(folder1+file, masked=False)
    print(src.rio.bounds(), src.rio.width, src.rio.height)
    src_files_to_mosaic.append(src)

mosaic = merge_arrays(src_files_to_mosaic)

print("The CRS:", mosaic.rio.crs)
print('Resolution:', mosaic.rio.resolution())
print("Bounds:", mosaic.rio.bounds())
print('Width:', mosaic.rio.width)
print('Height:', mosaic.rio.height)

(6.99986111111111, 44.999861111111116, 8.000138888888888, 46.0001388888889) 3601 3601
(7.99986111111111, 44.999861111111116, 9.00013888888889, 46.0001388888889) 3601 3601
(8.99986111111111, 44.999861111111116, 10.00013888888889, 46.0001388888889) 3601 3601
(6.99986111111111, 45.999861111111116, 8.000138888888888, 47.0001388888889) 3601 3601
(7.99986111111111, 45.999861111111116, 9.00013888888889, 47.0001388888889) 3601 3601
(8.99986111111111, 45.999861111111116, 10.00013888888889, 47.0001388888889) 3601 3601

The CRS: EPSG:4326
Resolution: (0.00027777777777777805, -0.0002777777777777778)
Bounds: (6.99986111111111, 44.999861111111116, 10.000138888888891, 47.0001388888889)
Width: 10801
Height: 7201

folder = 'folder2/'
files = ['ASTGTMV003_N45E007_dem_mask.tif', 'ASTGTMV003_N45E008_dem_mask.tif', 'ASTGTMV003_N45E009_dem_mask.tif',
        'ASTGTMV003_N46E007_dem_mask.tif', 'ASTGTMV003_N46E008_dem_mask.tif', 'ASTGTMV003_N46E009_dem_mask.tif']

src_files_to_mosaic = []
for file in files:
    src = rioxarray.open_rasterio(folder2+file, masked=False)
    print(src.rio.bounds(), src.rio.width, src.rio.height)
    src_files_to_mosaic.append(src)

mosaic_mask = merge_arrays(src_files_to_mosaic) 

print("The CRS:", mosaic_mask.rio.crs)
print('Resolution:', mosaic_mask.rio.resolution())
print("Bounds:", mosaic_mask.rio.bounds())
print('Width:', mosaic_mask.rio.width)
print('Height:', mosaic_mask.rio.height)

(6.99986111111111, 44.999861111111116, 8.000138888888888, 46.0001388888889) 3601 3601
(7.99986111111111, 44.999861111111116, 9.00013888888889, 46.0001388888889) 3601 3601
(8.99986111111111, 44.999861111111116, 10.00013888888889, 46.0001388888889) 3601 3601
(6.99986111111111, 45.999861111111116, 8.000138888888888, 47.0001388888889) 3601 3601
(7.99986111111111, 45.999861111111116, 9.00013888888889, 47.0001388888889) 3601 3601
(8.99986111111111, 45.999861111111116, 10.00013888888889, 47.0001388888889) 3601 3601

The CRS: EPSG:4326
Resolution: (0.00027777777777777805, -0.0002777777777777778)
Bounds: (6.99986111111111, 44.99958333333334, 10.000416666666668, 47.0001388888889) # <-- why these bounds not as before ?
Width: 10802 # <-- I expect this to be 10801
Height: 7202 # <-- I expect this to be 7201

The only way to have mosaic_mask have the expected shape 10801*7201 would be to pass the bounds argument as that of the previous mosaic. The documentation says that "If not set, bounds are determined from bounds of input DataArrays.". These bounds are the same in the two cases so I don't understand why if I don't specify them the resulting mosaic_mask has one extra col and one extra row.

mosaic_mask = merge_arrays(src_files_to_mosaic, bounds=mosaic.rio.bounds()) # <-- this would correctly yield a 10801 x 7201
print("Bounds:", mosaic_mask.rio.bounds())
print('Width:', mosaic_mask.rio.width)
print('Height:', mosaic_mask.rio.height)
Bounds: (6.99986111111111, 44.999861111111116, 10.000138888888891, 47.0001388888889)
Width: 10801
Height: 7201

[snowman2]

There are likely minor differences in the decimal precision of the transform/boundary your files that could cause the issue. I think the solution you provided is a great way to go.

[nmaffe]

The strage thing is that if the same merging exercise is done using rasterio instead of rioxarray the resulting mosaic_mask produced with rasterio is 10801 x 7201, so I'm not sure that the reason is the decimals.

[snowman2]

Mind printing the transform similar to how you did the bounds/resolution?

[nmaffe]

Thanks for the help @snowman2 !!

print("Bounds:", mosaic.rio.bounds())
print('Width:', mosaic.rio.width)
print('Height:', mosaic.rio.height)
print('Transform:', mosaic.rio.transform)

Bounds: (6.99986111111111, 44.999861111111116, 10.000138888888891, 47.0001388888889) Width: 10801 Height: 7201 Transform: <bound method XRasterBase.transform of <rioxarray.raster_array.RasterArray object at 0x7f6adc498910>>

print("Bounds:", mosaic_mask.rio.bounds())
print('Width:', mosaic_mask.rio.width)
print('Height:', mosaic_mask.rio.height)
print('Transform:', mosaic_mask.rio.transform)

Bounds: (6.99986111111111, 44.99958333333334, 10.000416666666668, 47.0001388888889) Width: 10802 Height: 7202 Transform: <bound method XRasterBase.transform of <rioxarray.raster_array.RasterArray object at 0x7f6ae722c070>>

I inspected the mosaic_mask (which consists of mostly zeros) and found out that the extra dimensions are one row to the bottom and one column to the right and both fully consist of ones (and this should not happens since the relevant tiles I'm merging do not contain ones the bottom and right borders), i.e.:

print(mosaic_mask)

<xarray.DataArray (band: 1, y: 7202, x: 10802)> array([[[0, 0, 0, ..., 0, 0, 1], [0, 0, 0, ..., 0, 0, 1], [0, 0, 0, ..., 0, 0, 1], ..., [0, 0, 0, ..., 0, 0, 1], [0, 0, 0, ..., 0, 0, 1], [1, 1, 1, ..., 1, 1, 1]]], dtype=int16)

So the problem comes down to why this ones are here.

[snowman2]

Thanks - as a side note, in rioxarray, transform is a method as there are options you can pass in. You will be able to see the value by calling it like: rio.transform(). Mind updating your response? Also, mind providing the transform output for the individual rasters?

[nmaffe]

Yeah, sure.

for file in files:
    src = rioxarray.open_rasterio(folder+file, masked=False)
    print(src.rio.bounds(), src.rio.width, src.rio.height)
    print(src.rio.transform())
    src_files_to_mosaic.append(src)
mosaic = merge_arrays(src_files_to_mosaic)
print("Bounds:", mosaic.rio.bounds())
print('Width:', mosaic.rio.width)
print('Height:', mosaic.rio.height)
print('Transform:', mosaic.rio.transform())

(6.99986111111111, 44.999861111111116, 8.000138888888888, 46.0001388888889) 3601 3601 | 0.00, 0.00, 7.00| | 0.00,-0.00, 46.00| | 0.00, 0.00, 1.00| (7.99986111111111, 44.999861111111116, 9.00013888888889, 46.0001388888889) 3601 3601 | 0.00, 0.00, 8.00| | 0.00,-0.00, 46.00| | 0.00, 0.00, 1.00| (8.99986111111111, 44.999861111111116, 10.00013888888889, 46.0001388888889) 3601 3601 | 0.00, 0.00, 9.00| | 0.00,-0.00, 46.00| | 0.00, 0.00, 1.00| (6.99986111111111, 45.999861111111116, 8.000138888888888, 47.0001388888889) 3601 3601 | 0.00, 0.00, 7.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00| (7.99986111111111, 45.999861111111116, 9.00013888888889, 47.0001388888889) 3601 3601 | 0.00, 0.00, 8.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00| (8.99986111111111, 45.999861111111116, 10.00013888888889, 47.0001388888889) 3601 3601 | 0.00, 0.00, 9.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00|

Bounds: (6.99986111111111, 44.999861111111116, 10.000138888888891, 47.0001388888889) Width: 10801 Height: 7201 Transform: | 0.00, 0.00, 7.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00|

for file in files_masked:
    src = rioxarray.open_rasterio(folder+file, masked=False)
    print(src.rio.bounds(), src.rio.width, src.rio.height)
    print(src.rio.transform())
    src_files_to_mosaic.append(src)    
mosaic_mask = merge_arrays(src_files_to_mosaic)
print("Bounds:", mosaic_mask.rio.bounds())
print('Width:', mosaic_mask.rio.width)
print('Height:', mosaic_mask.rio.height)
print('Transform:', mosaic_mask.rio.transform())

(6.99986111111111, 44.999861111111116, 8.000138888888888, 46.0001388888889) 3601 3601 | 0.00, 0.00, 7.00| | 0.00,-0.00, 46.00| | 0.00, 0.00, 1.00| (7.99986111111111, 44.999861111111116, 9.00013888888889, 46.0001388888889) 3601 3601 | 0.00, 0.00, 8.00| | 0.00,-0.00, 46.00| | 0.00, 0.00, 1.00| (8.99986111111111, 44.999861111111116, 10.00013888888889, 46.0001388888889) 3601 3601 | 0.00, 0.00, 9.00| | 0.00,-0.00, 46.00| | 0.00, 0.00, 1.00| (6.99986111111111, 45.999861111111116, 8.000138888888888, 47.0001388888889) 3601 3601 | 0.00, 0.00, 7.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00| (7.99986111111111, 45.999861111111116, 9.00013888888889, 47.0001388888889) 3601 3601 | 0.00, 0.00, 8.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00| (8.99986111111111, 45.999861111111116, 10.00013888888889, 47.0001388888889) 3601 3601 | 0.00, 0.00, 9.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00|

Bounds: (6.99986111111111, 44.99958333333334, 10.000416666666668, 47.0001388888889) Width: 10802 Height: 7202 Transform: | 0.00, 0.00, 7.00| | 0.00,-0.00, 47.00| | 0.00, 0.00, 1.00|

[snowman2]

That is definitely strange. Based on what you have posted, I don't see anything different between the rasters.