Merge Overlapping Rasters Using python and rioxarray

When utilizing tiled spatial data, it's common to come across overlapping tiles. For instance, when we chose four tiles from the global 30-meter land-cover dataset GLC_FCS30-2020 to encompass the entire North Island in New Zealand, there were overlapping pixels between them (Figure 1 and Figure 2).

The global 30-m land-cover data of GLC_FCS30-2020 are grouped by 948 5X5 degrees regional tiles in the GeoTIFF format, which are named “GLCFCS30_E/W**N/S**.tif”, where “E/W**N/S**” explains the longitude and latitude information of the upper left corner of each regional land-cover map. Further, each tile contains a land-cover label band ranging from 0–255, and the invalid fill value is labelled as 0 and 250.

Figure 1: Extracted 4 tiles for North Island in New Zealand with overlaps.

Figure 2: the right-bottom tile is dropped to show the overlap areas.

In the preceding tutorial, we utilized GDAL VRT Pixel Functions to merge overlapping tiles or rasters. In this tutorial, we aim to address the problem using another Python package called rioxarray, which expands upon xarray's geospatial capabilities with support from rasterio. Furthermore, we calculate the average value of overlapping pixels at the same location, rather than solely selecting maximum or minimum values from them as we did previously.

Solutions

import numpy as np
import xarray as xr
import rioxarray as rio
import matplotlib.pyplot as plt 

import warnings
warnings.filterwarnings('ignore')        

• Read tiles

tifs = ['data/GLC_FCS30_2020_E170S35.tif',     'data/GLC_FCS30_2020_E170S40.tif',
"data/GLC_FCS30_2020_E175S35.tif", 
"data/GLC_FCS30_2020_E175S40.tif"]

ls_data = []
for tif in tifs:
    ls_data.append(rio.open_rasterio(tif, mask_and_scale=True))        

• Merge tilesThe average is calculated as: average = sum/count. The sum and count functions are available from the rioxarray.merge.merge_arrays(), which are defined from rasterio. See rasterio.merge.merge() for details

raster_sum   = rio.merge.merge_arrays(dataarrays=ls_data, method='sum')
raster_count = rio.merge.merge_arrays(dataarrays=ls_data, method='count')

raster_avg = raster_sum/raster_count        

• Save data

raster_avg.rio.to_raster("rx_merged_GLC_FCS30_ave.tiff")         

Results

After merging these tiles, we still observe some strips which results from the overlapping sections (Figure 3). We will employ these sections to validate the accuracy of our merging process.

We can use GIS software to check the merged raster. When clicking on the strip part, we can confirm that the pixel value from the merged raster is indeed the average of the underlying four overlapped tiles at the same location (Figure 4).

Figure 3: Merged tiles. The strip showed the overlapped pixels.

Figure 4: Validation of pixel values on merged tiles against the selected tiles at a specific site.

Summary and Discussions

When conducting spatial analysis, it is often necessary to use tiled data, and encountering overlapping tiles is a common issue. Pixel functions or rioxarray/rasterio can be applied to solve the issue.

References

https://docs.xarray.dev/en/stable/

https://corteva.github.io/rioxarray/stable/

https://rasterio.readthedocs.io/en/stable/

https://www.dhirubhai.net/pulse/merge-overlapping-rasters-using-python-gdal-vrt-pixel-chonghua-yin-8ri1c?trackingId=HnZjg5Ifl%2FT3%2FpYmRdKvCA%3D%3D&lipi=urn%3Ali%3Apage%3Ad_flagship3_profile_view_base_recent_activity_content_view%3Bh%2FtDk02tSymsZDksThNNgw%3D%3D

Liangyun, L., Xiao, Z., Xidong, C., Yuan, G., & Jun, M. (2020). GLC_FCS30-2020: Global Land Cover with Fine Classification System at 30m in 2020 (v1.2) [Data set]. Zenodo.?https://doi.org/10.5281/zenodo.4280923