6

I'm using this code, of @Loïc Dutrieux, to stack 8 Landsat bands. 

with rasterio.open(final_list[0]) as src0:
   meta = src0.meta

with rasterio.open('stack.tif', 'w', **meta) as dst:
for id, layer in enumerate(final_list):
    with rasterio.open(layer) as src1:

        #the .astype(rasterio.int16) forces dtype to int16
        dst.write_band(id + 1, src1.read(1).astype(rasterio.int16))

When I open the stack.tif with gdalinfo I'm not getting bands' description in the Metadata: 

...
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
...
Band 1 Block=3596x1 Type=Int16, ColorInterp=Gray
NoData Value=-9999
Band 2 Block=3596x1 Type=Int16, ColorInterp=Undefined
NoData Value=-9999
...

I would like something similar to this:

...
Metadata:
Band_1=Band 1 Reflectance
Band_2=Band 2 Reflectance
Band_3=Band 3 Reflectance
Band_4=Band 4 Reflectance
Band_5=Band 5 Reflectance
Band_6=Band 7 Reflectance
Band_7=Band 6 Temperature
Band_8=Fmask
Image Structure Metadata:
...
Band 1 Block=300x1 Type=Int16, ColorInterp=Undefined
NoData Value=-9999
Band 2 Block=300x1 Type=Int16, ColorInterp=Undefined
NoData Value=-9999
...
mlateb
  • 139
  • 1
  • 6

1 Answers1

10

If using rasterio >= 1.0, use the dataset.set_band_description(self, bidx, value) method and dataset.descriptions property.

Sets the description of a dataset band.
Parameters
----------
bidx : int
    Index of the band (starting with 1).
value: string
    A description of the band.

For example:

descriptions = [
    'Band 1 Reflectance', 
    'Band 2 Reflectance',
    'Band 3 Reflectance',
    'Band 4 Reflectance',
    'Band 5 Reflectance',
    'Band 7 Reflectance', 
    'Band 6 Temperature'
]
for id, layer in enumerate(final_list, start=1):
    with rasterio.open(layer) as src1:

        #the .astype(rasterio.int16) forces dtype to int16
        dst.write_band(id, src1.read(1).astype(rasterio.int16))

        # Assuming the source dataset band 1 has a description
        dst.set_band_description(id, src1.descriptions[0])

        # If not, you can put anything in you want
        dst.set_band_description(id, descriptions[id-1])
user2856
  • 65,736
  • 6
  • 115
  • 196
  • Thank you @Luke for your answer. My rasterio.__version__ = '0.36.0'. Properties like src.descriptions are not available in my version. if it's a matter of version, please help me to install the newest. Could find any guide in the web – mlateb May 25 '18 at 21:55
  • Yes, not available in 0.36. Rasterio 1.0a12 is available if using conda in the conda-forge dev channel - conda install -c conda-forge/label/dev rasterio. I can't help with any other installation as I only use conda. – user2856 May 25 '18 at 22:06
  • If you can't upgrade rasterio, you can set the band descriptions with GDAL, but that's another question. Also, look into the rasterio dataset tags property, but I don't have rasterio 0.36 to check. – user2856 May 25 '18 at 22:10
  • Hi @Luke, I installed and used successfully Rasterio 1.0b1. The dst.set_band_description(id+1 , src1.descriptions[0]) added the band description as following: Band 1 Block=3596x1 Type=Int16, ColorInterp=Gray Description = band 2 surface reflectance NoData Value=-9999. But couln't add info to metadata section. I guess it's due to GTtiff driver specification. I opened a new question about masking and staking at the same time, please have a look. Thanks in advance! – mlateb May 30 '18 at 15:48
  • @mlateb Also have a look at the update_tags method http://rasterio.readthedocs.io/en/latest/topics/tags.html – user2856 May 30 '18 at 20:50