[gdal-dev] Hacking the new Landsat pixel_qa help...

[Rutger]

Jonathan Greenberg-4 wrote
> Folks:
> 
> USGS, in their infinite wisdom, decided to embed what is effectively a
> 1-bit, 16-band byte-interleaved-by-pixel image into a 16-bit, single-band
> geotiff for the Landsat quality mask.  I could write some crazy
> translation
> program that "reverse engineers" the data, but I was wondering if there is
> some way to "hack" or translate the file by simply redefining the "header"
> (change the number type to 1-bit, number of layers to 16, interleave to
> BIP).  Any suggestions for how to do this?

Hey Jonathan,

I personally think storing it in a single band is not that bad, because you
can use bitwise operations to test many flags in a single operation, which
is very efficient.

The new Python VRT pixel functions, which are available in 2.2 should be
able to do this. I haven't had a change to test it. But creating a VRT
calling some sort of 'bitcheck' function can be specified with a different
bit to check for each band. 

Modifying the example from the VRT tutorial on gdal.org, it might look
something like this (completely untested):

<VRTDataset rasterXSize="20" rasterYSize="20">
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
    <PixelFunctionType>bitcheck</PixelFunctionType>
    <PixelFunctionLanguage>Python</PixelFunctionLanguage>
    <PixelFunctionArguments check_bit="5"/>
    <PixelFunctionCode>
    </PixelFunctionCode>
    <SimpleSource>
      <SourceFilename
relativeToVRT="1">landsat_bqa_band.tif</SourceFilename>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

For multiple output bands it would be better to store the function in a
separate file somewhere, instead of repeating the code for each derived
band.

Since Landsat also combines bits (for the low/med/high ones), the function
would need to be changed a little, this is just to get an idea.  

I'm not convinced this is more better/preferable than doing the bit-checking
at the point where you would do something else with the data. Although maybe
for 'direct' visualization in QGIS for example it might be nice. 

The new VRT pixel functions open up a lot of possibilities, especially in
combination with Numba, its exciting to see where this is heading.


Regards,
Rutger
 




--
View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-Hacking-the-new-Landsat-pixel-qa-help-tp5319426p5319456.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.