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?