I have used gdalwarp for extracting Sentinel 5p
L2__AER_AI data product to 14 tiff files in EPSG:4326 CRS.
These projected tiff files can be gotten from mundi_tiff - Google Drive
Then, I used
gdal_merge.py for merging them:
gdal_merge.py -a_nodata 9.96920999999999943e+36 -o merged.tiff file_1.tiff file_2.tiff ...
The output looks correct, but it has pixel overlapping at the borders (opened on QGIS and applied style with pseudo color palette), and I don’t know what is the best way to flatten these pixels properly.