I am trying to create a 4-band raster image using gdal_merge.py but only the first band is getting the correct no data values assigned, according to ArcGIS and PCI. Here is my command:
gdal_merge.py -co COMPRESS=LZW -co BIGTIFF=YES -ul_lr %s -ot Int16 -separate -n -999 -a_nodata -999 -o %s %s' % (extent, outputFile, inputFilesString)
inputFilesString = "input1.tif input2.tif input3.tif input4.tif"
Where there is NoData, ArcGIS is showing NoData for Band1, but 0's for Bands2-4. However, this problem doesn't occur in QGIS (where all bands show NoData where they should). Also, according to gdalinfo (and ArcGIS/QGIS properties), each band's NoData value is correctly -999. All these things were leading me to think it could just be a bug in ArcGIS/PCI. However, when I read these bands (2-4) into individual arrays and write each of them each back to their own single band .tif (not setting the ND val) using GDAL API, it is writing these NoData areas as 0 while for band1 it correctly writes them as -999. So it is likely an issue in GDAL, not Arc/PCI.
I double checked (in ArcGIS and QGIS), and each of the individual tiffs that are the inputs for gdal_merge have NoData values of -999 that are correctly in areas where there should be NoData.
This is just not translating to the tiff output from gdal_merge.
Is this a bug with GDAL? Or am I not correctly assigning the NoData values in the gdal_merge.py command? From the documentation, it appears the NoData values should be assigned once for all bands as opposed to explicitly defining the values for each band, but I could be wrong.
UPDATE: I reran the exact same command as above but without the "-n -999" string, as suggested in an answer below, but nothing changed.