How to convert a Rasterio merge command into a GDAL merge command

Question

My understanding is that Rasterio is mostly a python wrapper for GDAL functions to make them easier to use. That certainly seems the case for me because I'm having a lot of trouble duplicating with GDAL what I was able to do easily with Rasterio's merge function.

I have two maps of roughly the same region that partially overlap, I want to use the lon,lat bounds and resolution of the first map image, and draw a second map image over it in the appropriate place.

In Rasterio it worked like this:

overlayMap,_ = merge([areaData,focusMapData], bounds=areaData.bounds, method='last')

Unfortunately this gives me a out-of-memory error in a few cases. I've read that to avoid that you should use GDAL directly with a virtual thingy like this:

areaData = rasterio.open('areaMap'+'.tif')
H, W = areaData.shape  ##- the order is switched from normal
bds= [areaData.bounds[i] for i in range(len(areaData.bounds))]
vrt = gdal.BuildVRT("merged.vrt", ['areaMap'+'.tif', 'focusMap'+'.tif'])
overlayMapData = gdal.Translate('/vsimem/in_memory_output.tif', vrt, outputBounds=bds, width=W, height=H)
vrt = None
overlayMap = overlayMapData.ReadAsArray()

But this doesn't work.

I get an output of the right size, but the focusMap data has been stretched to fill the area of the larger areaMap instead of being drawn on top of the areaMap matching its resolution in the appropriate part of the areaMap (as it is with the rasterio command).

So, what am I doing wrong? How do I get GDAL in python to do what I could already do easily in Rasterio?

Answer 1

The answer required two changes: (1) fixing the order of the bounds to match GDAL's top-left, bottom-right order, and (2) using the projWin command instead of the outputBounds command.

areaData = rasterio.open('areaMap'+'.tif')
H,W = areaData.shape  #- numpy output order is switched from "normal"
bds= [areaData.bounds[i] for i in range(len(areaData.bounds))]
bds = [bds [0], bds [3], bds [2], bds [1]] #- switch to GDAL's order
vrt = gdal.BuildVRT("merged.vrt", ['areaMap'+'.tif', 'focusMap'+'.tif'])
overlayMapData = gdal.Translate('/vsimem/in_memory_output.tif', vrt, 
projWin=bds, width=W, height=H)
vrt = None
overlayMap = overlayMapData.ReadAsArray()
overlayMap = np.dstack([overlayMap[0], overlayMap[1], overlayMap[2]])

Now the final overlapMap is in an array form that can be used with openCV and various AI methods.