Random Forest Classification

Install dependencies

%%bash
conda activate grass
conda install scikit-learn pandas

Setup

# Import modules
import io
import sys
import subprocess

# Find GRASS Python packages
sys.path.append(
  subprocess.check_output(
    ["grass", "--config", "python_path"],
    text=True
    ).strip()
  )

# Import GRASS packages
import grass.script as gs
import grass.jupyter as gj
from grass.tools import Tools

# Initialize GRASS project
gj.init("hilltop")
tools = Tools(overwrite=True)

# Print projection
projection = tools.g_proj(format="shell", flags="p").keyval
print(projection.get("name"))
NAD83(HARN) / Louisiana South

Define categories

# Define landcover categories
categories = """\
11|Water
24|Building
31|Bare
42|Forest
52|Shrub/Meadow
71|Grass
"""

Set region

# Set region
tools.g_region(raster="imagery", res=0.3)

Install extension

tools.g_extension(extension="r.learn.ml2")

Digitize training data

# Digitize
tools.g_gui_rdigit(create="training", type="CELL")

# Set landcover categories
tools.r_category(map="training", separator="pipe", rules=io.StringIO(categories))

# Set landcover colors
tools.r_colors(map="training", color="nlcd")

# Visualize
m = gj.Map(width=400)
m.d_rast(map="training")
m.show()

Random Forest

Supervised classification and regression with scikit-learn

# Set number of estimators
n = 600

# Remove image group
tools.g_remove(type="group", name="imagery", flags="f")

# Create image group
tools.i_group(group="imagery", input=["imagery.1", "imagery.2", "imagery.4", "imagery.4"])
tools.i_group(group="imagery", input=["ndvi", "ndwi"])
tools.i_group(group="imagery", input=["canopy"])
# tools.i_group(group="imagery", input=["landcover"])

# Training
tools.r_learn_train(
    group="imagery",
    training_map="training",
    model_name="RandomForestClassifier",
    n_estimators=n,
    save_model="random-forest.gz"
    )

# Prediction
tools.r_learn_predict(
    group="imagery",
    load_model="random-forest.gz",
    output="random_forest"
    )

# Set landcover colors
tools.r_colors(map="random_forest", color="nlcd")

# Set landcover categories
tools.r_category(map="random_forest", separator="pipe", rules=io.StringIO(categories))

# Visualize landcover
m = gj.Map(width=800)
m.d_rast(map="random_forest")
m.d_legend(
    raster="random_forest",
    color="white",
    at=(75, 95, 1, 3),
    use=(11,24,31,42,52,71))
m.show()