2

I am trying the following code to mosaic/merge 10 image files (.tif files) using the following code:

import rioxarray
from rioxarray import merge
from rasterio.plot import show

items = [r'D:\LC08_L1TP_146046_20180425_20180502_01_T1.tif',
r'D:\LE07_L1TP_145045_20180426_20180523_01_T1.tif',
r'D:\LC08_L1TP_144046_20180427_20180502_01_T1.tif',
r'D:\LC08_L1TP_144045_20180427_20180502_01_T1.tif',
r'D:\LC08_L1TP_143046_20180420_20180502_01_T1.tif',
r'D:\LC08_L1TP_143047_20180420_20180502_01_T1.tif',
r'D:\LC08_L1TP_143045_20180420_20180502_01_T1.tif',
r'D:\LE07_L1TP_145046_20180426_20180523_01_T1.tif',
r'D:\LE07_L1TP_145045_20180426_20180523_01_T1.tif',
r'D:\LC08_L1TP_144047_20180427_20180502_01_T1.tif']

elements = []

for val in items:
    elements.append(rioxarray.open_rasterio(val))

merged = merge.merge_arrays(elements, nodata=0.0)

image = merged.values
show(image)

After running the above code I am getting the following plot: Image after running above code

The image shown above is not the correct mosaic/merge image as the left and right part of the image is swapped because when I am loading the input files in ArcMap I am getting the following:

Input Files as seen in ArcMap

How can I resolve this issue?

2
  • 1
    It looks like two different UTM zones are spanned, so the western half is east of the eastern half.
    – Vince
    Oct 15, 2020 at 11:28
  • @Vince. Your suggestion was correct. It was because of Projection difference (difference in UTM zones between input images) because of which mosaic was not happening properly.
    – RRSC NGP
    Oct 16, 2020 at 4:21

1 Answer 1

3

Cross posted: https://github.com/corteva/rioxarray/issues/172

I see that you have both EPSG:32643 and EPSG:32644 as the CRS of the datasets. The merge function assumes that they are all the same CRS.

You need to change:

for val in items:
    elements.append(rioxarray.open_rasterio(val))

to:

for val in items:
    rds = rioxarray.open_rasterio(val)
    if rds.rio.crs != "EPSG:32643":
        rds = rds.rio.reproject("EPSG:32643")
    elements.append(rds)
2
  • After using your suggestion I got following error: "AttributeError: 'DataArray' object has no attribute 'crs'". So changing "rds.crs" to "rds.rio.crs" solve the issue. Thanks for your help.
    – RRSC NGP
    Oct 16, 2020 at 4:19
  • Good catch, glad you were able to find the fix 👍
    – snowman2
    Oct 16, 2020 at 12:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.