In this example we will demonstrate how to fit and score a supervised learning model with a sample of Sentinel-2 data and hand-drawn vector labels over different [land cover](https://en.wikipedia.org/wiki/Land_cover) types.
```python, setup, echo=False
from IPython.core.display import display
from pyrasterframes.utils import create_rf_spark_session
from pyrasterframes.rasterfunctions import *
import pyrasterframes.rf_ipython
from pyspark.sql.functions import lit
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
spark = create_rf_spark_session()
```
## Create and Read Raster Catalog
The first step is to create a Spark DataFrame containing our imagery data. To achieve that we will create @ref:[a catalog DataFrame](raster-catalogs.md#creating-a-catalog). In the catalog, each row represents a distinct area and time, and each column is the URI to a band's image product. In this example our catalog just has one row. After reading the catalog, the resulting Spark DataFrame may have many rows per URI, with a column corresponding to each band.
The imagery for feature data will come from [eleven bands of 60 meter resolution Sentinel-2](https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/resolutions/spatial) imagery. We also will use the [scene classification (SCL)](https://earth.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm) data to identify high quality, non-cloudy pixels.
The land classification labels are based on a small set of hand drawn polygons in the GeoJSON file [here](https://github.com/locationtech/rasterframes/blob/develop/pyrasterframes/src/test/resources/luray-labels.geojson). The property `id` indicates the type of land cover in each area. For these integer values, 1 is forest, 2 is cropland, and 3 is developed areas.
We will create a very small Spark DataFrame of the label shapes and then join it to the raster DataFrame. Such joins are typically expensive, but in this case both datasets are quite small. To speed up the join for the small vector DataFrame, we put the `broadcast` hint on it, which will tell Spark to put a copy of it on each Spark executor.
After the raster and vector data are joined, we will convert the vector shapes into _tiles_ using the @ref:[`rf_rasterize`](reference.md#rf-rasterize) function. This procedure is sometimes called "burning in" a geometry into a raster. The values in the resulting _tile_ cells are the `id` property of the GeoJSON, which we will use as labels in our supervised learning task. In areas where the geometry does not intersect, the cells will contain NoData.
To filter only for good quality pixels, we follow roughly the same procedure as demonstrated in the @ref:[quality masking](nodata-handling.md#masking) section of the chapter on NoData. Instead of actually setting NoData values in the unwanted cells of any of the imagery bands, we will just filter out the mask cell values later in the process.
We import various Spark components that we need to construct our [`Pipeline`](https://spark.apache.org/docs/latest/ml-pipeline.html). These are the objects that will work in sequence to conduct the data preparation and modeling.
```python, imports, echo=True
from pyrasterframes import TileExploder
from pyrasterframes.rf_types import NoDataFilter
from pyspark.ml.feature import VectorAssembler
from pyspark.ml.classification import DecisionTreeClassifier
from pyspark.ml.evaluation import MulticlassClassificationEvaluator
from pyspark.ml import Pipeline
```
SparkML requires that each observation be in its own row, and those observations be packed into a single [`Vector`](https://spark.apache.org/docs/latest/api/python/pyspark.ml.html#module-pyspark.ml.linalg) object. The first step is to "explode" the _tiles_ into a single row per cell or pixel with the `TileExploder` (see also @ref:[`rf_explode_tiles`](reference.md#rf_explode_tiles)). If a _tile_ cell contains a NoData it will become a null value after the exploder stage. Then we use the `NoDataFilter` to filter out any rows that missing or null values, which will cause an error during training. Finally we use the SparkML `VectorAssembler` to create that `Vector`.
Recall above we set undesirable pixels to NoData, so the `NoDataFilter` will remove them at this stage. We apply the filter to the `mask` column and the `label` column, the latter being used during training. When it is time to score the model, the pipeline will ignore the fact that there is no `label` column on the input DataFrame.
```python, transformers
exploder = TileExploder()
noDataFilter = NoDataFilter() \
.setInputCols(['label', 'mask'])
assembler = VectorAssembler() \
.setInputCols(bands) \
.setOutputCol("features")
```
We are going to use a decision tree for classification. You can swap out one of the other multi-class classification algorithms if you like. With the algorithm selected we can assemble our modeling pipeline.
The next step is to actually run each step of the Pipeline we created, including fitting the decision tree model. We filter the DataFrame for only _tiles_ intersecting the label raster because the label shapes are relatively sparse over the imagery. It would be logically equivalent to either include or exclude thi step, but it is more efficient to filter because it will mean less data going into the pipeline.
To view the model's performance, we first call the pipeline's `transform` method on the training dataset. This transformed dataset will have the model's prediction included in each row. We next construct an evaluator and pass it the transformed dataset to easily compute the performance metric. We can also create custom metrics using a variety of DataFrame or SQL transformations.
```python, eval
prediction_df = model.transform(df_mask) \
.drop(assembler.getOutputCol()).cache()
prediction_df.printSchema()
eval = MulticlassClassificationEvaluator(
predictionCol=classifier.getPredictionCol(),
labelCol=classifier.getLabelCol(),
metricName='accuracy'
)
accuracy = eval.evaluate(prediction_df)
print("\nAccuracy:", accuracy)
```
As an example of using the flexibility provided by DataFrames, the code below computes and displays the confusion matrix.
Because the pipeline included a `TileExploder`, we will recreate the tiled data structure. The explosion transformation includes metadata enabling us to recreate the _tiles_. See the @ref:[`rf_assemble_tile`](reference.md#rf-assemble-tile) function documentation for more details. In this case, the pipeline is scoring on all areas, regardless of whether they intersect the label polygons. This is simply done by removing the `label` column, as @ref:[discussed above](supervised-learning.md#create-ml-pipeline).
Take a look at a sample of the resulting prediction and the corresponding area's red-green-blue composite image. Note that because each `prediction` tile is rendered independently, the colors may not have the same meaning across rows.