Creating multi-image single band single file geotiff using GDAL with NumPy?

Question

I have a geotiff file that is single band but contains multiple subdatasets. I use gdal with numpy to open a specific subdataset and mask some of its values.

Example:

$ gdalinfo sample.tif
...
Subdatasets:
SUBDATASET_1_NAME=GTIFF_DIR:1:sample.tif
SUBDATASET_1_DESC=...
SUBDATASET_2_NAME=GTIFF_DIR:2:sample.tif
SUBDATASET_2_DESC=...

# pseudocode
sds = gdal.Open('GTIFF_DIR:2:sample.tif')
sds_mask = apply mask to sds

I want to put sds_mask back into the original geotiff or a copy of it. Is this possible with gdal?

I tried writing out all the subdatasets, modifying one of the subdatasets, and then merging them back together with gdal_translate and gdal_merge.py with the -separate option but it does not create a geotiff with subdatasets.

Answer 1

GDAL can't do this to my knowledge, although LibTIFF (http://www.simplesystems.org/libtiff/) has a tool called "tiffcp" that can do this:

To create a file out.tif that has two subdatasets a.tif and b.tif:

tiffcp a.tif b.tif out.tif

If the TIF will be larger than 4GB, you'll have to use the "-8" option