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

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Jonathan Greenberg-4
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?  

I dropped an example file here:


To be clear, the "displayed" number has no real meaning, the information is embedded into the set of individual bits in each of the 16 positions.

--j
--
-- 
Jonathan A. Greenberg, PhD
Randall Endowed Professor and Associate Professor of Remote Sensing
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Natural Resources & Environmental Science
University of Nevada, Reno
1664 N Virginia St MS/0186
Reno, NV 89557
Phone: 415-763-5476
http://www.unr.edu/nres
AIM: jgrn307, MSN: [hidden email], Gchat: jgrn307, Skype: jgrn3007

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: [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>
<![CDATA[
        import numpy as np
        def bitcheck(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                     raster_ysize, buf_radius, gt, **kwargs):
        
            check_bit = int(kwargs['check_bit'])
            out_ar = (in_ar & check_bit) > 0
    ]]>

    </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
 
SBL
Reply | Threaded
Open this post in threaded view
|

Re: Hacking the new Landsat pixel_qa help...

SBL
Hi,

If you don't mind some advertising of own work, have a look at:
https://grass.osgeo.org/grass72/manuals/addons/i.landsat8.qc.html

If you want to do this using gdal,  the python code might be a source of inspiration...

Cheers,
Stefan
________________________________________
Von: gdal-dev [[hidden email]] im Auftrag von Rutger [[hidden email]]
Gesendet: Freitag, 5. Mai 2017 10:37
An: [hidden email]
Betreff: Re: [gdal-dev] Hacking the new Landsat pixel_qa help...

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.
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Hacking the new Landsat pixel_qa help...

Matthew Perry-2
In reply to this post by Jonathan Greenberg-4
Hi Jonathan,

Jonathan Greenberg-4 wrote
> 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?


https://github.com/mapbox/rio-l8qa might be what you're looking for; provides a python API and CLI for extracting arrays out of the Landsat QA band, both the pre-Collection and collection format. It's at version 0.1 so caveat emptor but could be useful. Let me know how it works for you.

Matt Perry


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev