[gdal-dev] keep metadata tags when creating multiband file

[Even Rouault]

Michele,

> 
> I have a number of individual Landsat bands that I'm converting to GeoTiff
> and would then like to merge to a multiband file using either
> gdal_merge.py or another gdal utility (gdalbuildvrt, gdal_translate).  In
> the individually derived geotif files, I am able to internally "label" the
> corresponding Landsat band by creating a metadata tag.  For example, in
> using the following gdal_translate command from a file that is Landsat
> data for Band 1 called landsat.b1.dat:
> 
> gdal_translate -of GTiff -strict -stats -mo "Band = LandsatBand ${f: -5:1}"
> -co "COMPRESS=LZW"  -co "PROFILE=GeoTIFF" $f $f.tif
> 
> This creates a metadata tag in the tif file:
> -------------------------------------------------
> Metadata:
>   AREA_OR_POINT=Area
>   Band =LandsatBand 1
> 
> And there is a separate .aux.xml for each of the tif files:
> -------------------------------------------------------------------
> <PAMDataset>
>   <Metadata>
>     <MDI key="Band ">LandsatBand 1</MDI>
>   </Metadata>
>   <PAMRasterBand band="1">
>     <Metadata>
>       <MDI key="STATISTICS_MAXIMUM">255</MDI>
>       <MDI key="STATISTICS_MEAN">40.826976974681</MDI>
>       <MDI key="STATISTICS_MINIMUM">0</MDI>
>       <MDI key="STATISTICS_STDDEV">31.471508896152</MDI>
>     </Metadata>
>   </PAMRasterBand>
> 

Instead of setting the band name at the dataset level, you should set it at 
the band level instead. There's currently no way to do that with 
gdal_translate. You have to edit the .aux.xml file manually :

For example let's consider band1.tif with band1.tif.aux.xml :

<PAMDataset>
  <PAMRasterBand band="1">
    <Metadata>
      <MDI key="Band">LandsatBand 1</MDI>
    </Metadata>
  </PAMRasterBand>
</PAMDataset>

and band2.tif with band2.tif.aux.xml :

<PAMDataset>
  <PAMRasterBand band="1">
    <Metadata>
      <MDI key="Band">LandsatBand 2</MDI>
    </Metadata>
  </PAMRasterBand>
</PAMDataset>

and then you would want to run "gdalbuildvrt -separate merged.vrt band1.tif 
band2.tif". Unfortunately I see that it also loses band metadata. I don't 
think there's a fundamental reason to explain that they are lost, except that 
it hasn't been coded yet.

You can still put then again by editing the merged.vrt and adding the     
<Metadata>      <MDI key="Band">LandsatBand X</MDI>     </Metadata> at 
appropriate place.

See below :

<VRTDataset rasterXSize="20" rasterYSize="20">
  <SRS>PROJCS["NAD27 / UTM zone 
11N",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 
1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","26711"]]</SRS>
  <GeoTransform>  4.4072000000000000e+05,  6.0000000000000000e+01,  
0.0000000000000000e+00,  3.7513200000000000e+06,  0.0000000000000000e+00, 
-6.0000000000000000e+01</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1">
    <Metadata>
      <MDI key="Band">LandsatBand 1</MDI>
    </Metadata>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">byte1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="20" RasterYSize="10" DataType="Byte" 
BlockXSize="20" BlockYSize="10" />
      <SrcRect xOff="0" yOff="0" xSize="20" ySize="10" />
      <DstRect xOff="0" yOff="0" xSize="20" ySize="10" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="2">
    <Metadata>
      <MDI key="Band">LandsatBand 2</MDI>
    </Metadata>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">byte2.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="20" RasterYSize="20" DataType="Byte" 
BlockXSize="20" BlockYSize="20" />
      <SrcRect xOff="0" yOff="0" xSize="20" ySize="20" />
      <DstRect xOff="0" yOff="0" xSize="20" ySize="20" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>


And then

$ gdal_translate merged.vrt merged.tif 

$ gdalinfo merged.tif 
Driver: VRT/Virtual Raster
Files: merged.tif 
Size is 20, 20
Coordinate System is:
PROJCS["NAD27 / UTM zone 11N",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982139006,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26711"]]
Origin = (440720.000000000000000,3751320.000000000000000)
Pixel Size = (60.000000000000000,-60.000000000000000)
Corner Coordinates:
Upper Left  (  440720.000, 3751320.000) (117d38'28.21"W, 33d54' 8.47"N)
Lower Left  (  440720.000, 3750120.000) (117d38'27.92"W, 33d53'29.51"N)
Upper Right (  441920.000, 3751320.000) (117d37'41.48"W, 33d54' 8.71"N)
Lower Right (  441920.000, 3750120.000) (117d37'41.20"W, 33d53'29.75"N)
Center      (  441320.000, 3750720.000) (117d38' 4.70"W, 33d53'49.11"N)
Band 1 Block=20x20 Type=Byte, ColorInterp=Gray
  Metadata:
    Band=LandsatBand 1
Band 2 Block=20x20 Type=Byte, ColorInterp=Undefined
  Metadata:
    Band=LandsatBand 2

Best rergards,

Even

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html