[GRASS-user] help to use r.mapcalc with decimal/binary (modisdata)

[Glynn Clements]

andrew.haywood at poyry.com wrote:

> Okay - unfortunately i am now totally confused. I have tried to read up on 
> binary data (to little or no avail!!!). My understanding is the MODIS 
> suf_refl_qc data that i have import is a unsigned 16 bit integer number. I 
> would like to select out pixels using this raster  to create a cloud mask 
> based on some of the bit codes.
> 
> for example identify " band 1 quality = missing input" would have the 
> following 
> 
>  * * * * 1 0 1 1 * * * * * * * * * - where * represents placeholder
> 
> If i try and follow Glynn's logic he is saying divide the number by 16 and 
> then take the modulus 16 and then equate it to "11" to get the subset. 
> 
> I am unsure how I choose the "11" and how do i change the syntax to pull 
> out the different masks i may be interested in.
> 
> 0 0 * * * * * * * * * * * * * *   ideal quality all bands
> 0 1 * * * * * * * * * * * * * * less than ideal quality
> * * 0 0 * * * * * * * * * * * * clear
> * * 0 1 * * * * * * * * * * * * cloudy
> * * 1 0 * * * * * * * * * * * * mixed
> * * * * 0 0 0 0* * * * * * * * highest quality band 1
> * * * * * * * * * 0 0 0 0* * * highest quality band 2
> etc...
> 
> Any help would be greatly appreciated. At this stage im plugging in the 
> decimal numbers into binary calculator and using r.reclass to create a 
> MASK.

First, binary numbers are normally written with the bit 0 (the "units"
digit) on the right, so "band 1 quality = missing input" would be:

	1 1 1 1 1 1
	5 4 3 2 1 0 9 8 7 6 5 4 3 2 1 0

	* * * * * * * * 1 0 1 1 * * * *

To shift everything right by N bits, divide by 2^N, so dividing by 16
(2^4) shifts everything right by 4 places, moving bits 4-7 to bits
0-3:

	* * * * * * * * * * * * 1 0 1 1

[Similarly, dividing a decimal number by 10^4 = 10,000 shifts it right
by 4 digits: 12345678 / 10000 = 1234]

To keep the bottom N bits and discard the rest, take the remainder
(modulus) after division by 2^N, so modulus 16 (2^4) keeps the bottom
4 bits:

	0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1

[Similarly dividing a decimal number by 10^4 = 10,000 and taking the
remainder takes the bottom 4 digits: 12345678/10000 = 1234 remainder
5678, so 12345678 % 10000 = 5678]

The binary number "1011"
	= 1*2^3 + 0*2"2 + 1*2^1 + 1*2^0
	= 1*8 + 0*4 + 1*2 + 1*1
	= 8 + 2 + 1
	= 11

So, given your original bit definitions

	bits	description		divide	modulus	result range
	0-1	MODLAND QA		1	4	0-3
	2-3	cloud state		4	4	0-3
	4-7	band 1 data quality	16	16	0-15
	8-11	band 2 data quality	256	16	0-15
	12	atmospheric correction	4096	2	0-1
	13	adjacency correction	8192	2	0-1
	14	different orbit		16384	2	0-1
	15	unused			32768	2	0-1

For the right-hand side:

MODLAND QA (2 bits):

	binary	decimal	description
	00	 0	corrected product produced at ideal quality -- all bands
	01	 1	corrected product produced, less than ideal quality -- some or all bands
	10	 2	corrected product not produced due to cloud effects -- all bands
	11	 3	corrected product not produced for other reasons -- some or all

Cloud state (2 bits):

	binary	decimal	description
	00	 0	clear		    
	01	 1	cloudy		    
	10	 2	mixed		    
	11	 3	not set, assumed clear

band quality (4 bits):

	binary	decimal	description
	0000	 0	highest quality
	0001	 1
	0010	 2
	0011	 3
	0100	 4
	0101	 5
	0110	 6
	0111	 7
	1000	 8	dead detector; data has been copied from adjacent detector
	1001	 9	solar zenith >= 86 degrees										
	1010	10	solar zenith >= 85 and < 86 degrees									
	1011	11	missing input												
	1100	12	internal constant used in place of climatological data for at least one atmospheric constant		
	1101	13	correction out of bounds, pixel constrained to extreme allowable value					
	1110	14	L1B data faulty												
	1111	15	not processed due to deep ocean or clouds								

A shell script which uses r.mapcalc to split a QC map given as its
first argument into separate maps for each field:

#!/bin/sh
src=$1
r.mapcalc <<EOF
$src.modland    = ($src / 1    ) %  4
$src.cloud      = ($src / 4    ) %  4
$src.band1qual  = ($src / 16   ) % 16
$src.band2qual  = ($src / 256  ) % 16
$src.atmos_corr = ($src / 4096 ) %  2
$src.adj_corr   = ($src / 8192 ) %  2
$src.diff_orbit = ($src / 16384) %  2
$src.unused     = ($src / 32768) %  2
EOF

-- 
Glynn Clements <glynn at gclements.plus.com>