Add landsat masking section to masking docs page

Signed-off-by: Jason T. Brown <[email protected]>

Author: vpipkt

pyrasterframes/src/main/python/docs/masking.pymd

@@ -14,8 +14,10 @@ spark = create_rf_spark_session()
 ```
 
 Masking is a common operation in raster processing. It is setting certain cells to the @ref:[NoData value](nodata-handling.md). This is usually done to remove low-quality observations from the raster processing. Another related use case is to @ref:["clip"](masking.md#clipping) a raster to a given polygon.  
+
+In this section we will demonstrate two common schemes for masking. In Sentinel 2, there is a separate classification raster that defines low quality areas. In Landsat 8, several quality factors are measured and the indications are packed into a single integer, which we have to unpack.
 
-## Masking Example
+## Masking Sentinel 2
 
 Let's demonstrate masking with a pair of bands of Sentinel-2 data. The measurement bands we will use, blue and green, have no defined NoData. They share quality information from a separate file called the scene classification (SCL), which delineates areas of missing data and probable clouds. For more information on this, see the [Sentinel-2 algorithm overview](https://earth.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm). Figure 3 tells us how to interpret the scene classification. For this example, we will exclude NoData, defective pixels, probable clouds, and cirrus clouds: values 0, 1, 8, 9, and 10.
 
@@ -40,7 +42,7 @@ unmasked.printSchema()
 unmasked.select(rf_cell_type('blue'), rf_cell_type('scl')).distinct()
 ```
 
-## Define CellType for Masked Tile
+### Define CellType for Masked Tile
 
 Because there is not a NoData already defined for the blue band, we must choose one. In this particular example, the minimum value is greater than zero, so we can use 0 as the NoData value. We will construct a new `CellType` object to represent this. 
 
@@ -59,7 +61,7 @@ We next convert the blue band to this cell type.
 converted = unmasked.select('scl', 'green', rf_convert_cell_type('blue', masked_blue_ct).alias('blue'))
 ```
 
-## Apply Mask from Quality Band
+### Apply Mask from Quality Band
 
 Now we set cells of our `blue` column to NoData for all locations where the `scl` tile is in our set of undesirable values. This is the actual _masking_ operation.
 
@@ -89,7 +91,7 @@ And the original SCL data. The bright yellow is a cloudy region in the original
 display(sample[1])
 ```
 
-## Transferring Mask
+### Transferring Mask
 
 We can now apply the same mask from the blue column to the green column. Note here we have supressed the step of explicitly checking what a "safe" NoData value for the green band should be.  
 
@@ -98,9 +100,72 @@ masked.withColumn('green_masked', rf_mask(rf_convert_cell_type('green', masked_b
       .orderBy(-rf_no_data_cells('blue_masked'))
 ```
 
+## Masking Landsat 8
+
+
+We will work with the Landsat scene [here](https://landsat-pds.s3.us-west-2.amazonaws.com/c1/L8/153/075/LC08_L1TP_153075_20190718_20190731_01_T1/index.html). For simplicity, we will just use two of the seven 30m bands. The quality mask for all bands is all contained in the `BQA` band. 
+
+
+```python, build_l8_df 
+base_url = 'https://landsat-pds.s3.us-west-2.amazonaws.com/c1/L8/153/075/LC08_L1TP_153075_20190718_20190731_01_T1/LC08_L1TP_153075_20190718_20190731_01_T1_'
+data4 = base_url + 'B4.TIF'
+data2 = base_url + 'B2.TIF'
+mask = base_url + 'BQA.TIF'
+l8_df = spark.read.raster([[data4, data2, mask]]) \
+               .withColumnRenamed('proj_raster_0', 'data') \
+               .withColumnRenamed('proj_raster_1', 'data2') \
+               .withColumnRenamed('proj_raster_2', 'mask')
+```
+
+Masking is described [on the Landsat Missions page](https://www.usgs.gov/land-resources/nli/landsat/landsat-collection-1-level-1-quality-assessment-band). It is pretty dense. Focus for this data set is the Collection 1 Level-1 for Landsat 8. 
+ 
+There are several inter-related factors to consider. In this exercise we will mask away the following.
+ 
+ * Designated Fill = yes
+ * Cloud = yes
+ * Cloud Shadow Confidence = Medium or High
+ * Cirrus Confidence = Medium or High
+ 
+Note that you should consider your application and do your own exploratory analysis to determine the most appropriate mask!
+ 
+According to the information on the Landsat site this translates to masking by bit values in the BQA according to the following table.
+ 
+| Description        	| Value    	| Bits  	| Bit values     	|
+|--------------------	|----------	|-------	|----------------	|
+| Designated fill    	| yes      	| 0     	| 1              	|
+| Cloud              	| yes      	| 4     	| 1              	|
+| Cloud shadow conf. 	| med / hi 	| 7-8   	| 10, 11 (2, 3)  	|
+| Cirrus conf.       	| med / hi 	| 11-12 	| 10, 11 (2, 3)  	|
+
+
+In this case, we will use the value of 0 as the NoData in the band data. Inspecting the associated [MTL txt file](https://landsat-pds.s3.us-west-2.amazonaws.com/c1/L8/153/075/LC08_L1TP_153075_20190718_20190731_01_T1/LC08_L1TP_153075_20190718_20190731_01_T1_MTL.txt) By inspection, we can discover that the minimum value in the band will be 1, thus allowing our use of 0 for NoData.
+
+The code chunk below works through each of the rows in the table above. The first expression sets the cell type to have the selected NoData. The @ref:[`rf_mask_by_bit`](reference.md#rf-mask-by-bit) and @ref:[`rf_mask_by_bits`](reference.md#rf-mask-by-bits) functions extract the selected bit or bits from the `mask` cells and compare them to the provided values.
+
+```python, build_l8_mask
+l8_df = l8_df.withColumn('data_masked', # set to cell type that has a nodata 
+                   rf_convert_cell_type('data', CellType.uint16())) \
+       .withColumn('data_masked', # fill yes 
+                  rf_mask_by_bit('data_masked', 'mask', 0, 1)) \
+       .withColumn('data_masked', # cloud yes
+                  rf_mask_by_bit('data_masked', 'mask', 4, 1)) \
+       .withColumn('data_masked', # cloud shadow conf is medium or high
+                  rf_mask_by_bits('data_masked', 'mask', 7, 2, [2, 3])) \
+       .withColumn('data_masked', # cloud shadow conf is medium or high
+                  rf_mask_by_bits('data_masked', 'mask', 11, 2, [2, 3])) \
+       .withColumn('data2', # mask other data col against the other band
+                  rf_mask(rf_convert_cell_type('data2',CellType.uint16()), 'data_masked')) \
+       .filter(rf_data_cells('data_masked') > 0) # remove any entirely ND rows
+
+# Inspect a sample
+l8_df.select('data', 'mask', 'data_masked', 'data2', rf_extent('data_masked')) \
+     .filter(rf_data_cells('data_masked') > 32000)
+```
+
+
 ## Clipping
 
-Clipping is the use of a polygon to determine the areas to mask in a raster. Typically the areas inside a polygon are retained and the cells outside are set to NoData. Given a geometry column on our DataFrame, we have to carry out three basic steps. First we have to ensure the vector geometry is correctly projected to the same @ref:[CRS](concepts.md#coordinate-reference-system-crs) as the raster. We'll continue with our example creating a simple polygon. Buffering a point will create an approximate circle.
+Clipping is the use of a polygon to determine the areas to mask in a raster. Typically the areas inside a polygon are retained and the cells outside are set to NoData. Given a geometry column on our DataFrame, we have to carry out three basic steps. First we have to ensure the vector geometry is correctly projected to the same @ref:[CRS](concepts.md#coordinate-reference-system-crs) as the raster. We'll continue with our Sentinel 2 example, creating a simple polygon. Buffering a point will create an approximate circle.
 
 
 ```python, reproject_geom