Chapter 17: Advanced Spatial Modeling

Learning Objectives

After reading this chapter, you should be able to:

  • Understand advanced spatial constraints in forest optimization

  • Implement adjacency and contiguous area requirements

  • Model spatial connectivity and fragmentation

  • Use raster-based spatial analysis with ws3

  • Integrate GIS data into optimization models

Introduction

Spatial modeling is crucial for realistic forest management optimization. While basic ws3 models treat all areas as homogeneous, real forests have spatial structure that affects:

  • Harvest logistics: Distance to roads, landings, mills

  • Ecological requirements: Habitat connectivity, buffer zones

  • Regulatory constraints: Adjacency rules, contiguous areas

  • Visual impacts: Screen distances, viewshed analysis

This chapter explores advanced spatial modeling techniques that go beyond basic spatial allocation.

Spatial Constraints Overview

Common spatial constraints in forest optimization:

  1. Adjacency Constraints: Prevent harvesting adjacent areas

  2. Contiguous Area Requirements: Ensure harvested areas meet size minimums

  3. Spatial Connectivity: Maintain landscape connectivity

  4. Buffer Zones: Protect sensitive areas from harvesting

  5. Visual Quality: Manage visual impacts from harvest blocks

Mathematical Formulation:

For adjacency constraints, we typically use:

\[x_{i,t} + x_{j,t} \leq 1 \quad \forall i,j \in \text{adjacent pairs}, \forall t\]

Where \(x_{i,t}\) is a binary variable indicating whether area \(i\) is harvested in period \(t\).

Implementing Adjacency Constraints

Step 1: Define Adjacency Matrix

import numpy as np
import geopandas as gpd

# Load spatial data
spatial_df = gpd.read_file("data/spatial_inventory.geojson")

# Create adjacency matrix
def create_adjacency_matrix(geometry_series):
    """Create binary adjacency matrix from geometries."""
    n = len(geometry_series)
    adj_matrix = np.zeros((n, n), dtype=int)

    for i in range(n):
        for j in range(i+1, n):
            if geometry_series[i].touches(geometry_series[j]):
                adj_matrix[i, j] = 1
                adj_matrix[j, i] = 1

    return adj_matrix

adj_matrix = create_adjacency_matrix(spatial_df.geometry)
print(f"Adjacency matrix shape: {adj_matrix.shape}")
print(f"Number of adjacent pairs: {adj_matrix.sum() // 2}")

Step 2: Add Adjacency Constraints to Problem

from ws3.opt import Problem

def add_adjacency_constraints(problem, adj_matrix, dt_mapping):
    """Add adjacency constraints to optimization problem.

    :param problem: ws3 optimization problem
    :param adj_matrix: adjacency matrix
    :param dt_mapping: mapping from spatial units to DT codes
    """
    n_spatial = adj_matrix.shape[0]

    # Find adjacent pairs
    adjacent_pairs = []
    for i in range(n_spatial):
        for j in range(i+1, n_spatial):
            if adj_matrix[i, j] == 1:
                adjacent_pairs.append((i, j))

    print(f"Adding {len(adjacent_pairs)} adjacency constraints")

    # Add constraints
    for idx, (i, j) in enumerate(adjacent_pairs):
        # Get DT codes for adjacent areas
        dt_i = dt_mapping[i]
        dt_j = dt_mapping[j]

        # Create constraint: x_{i,t} + x_{j,t} <= 1
        constraint_name = f"adj_{i}_{j}"
        coeffs = {
            f"{dt_i}_t": 1.0,
            f"{dt_j}_t": 1.0,
        }

        problem.add_constraint(
            name=constraint_name,
            coeffs=coeffs,
            sense='leq',
            rhs=1.0
        )

Step 3: Solve with Adjacency Constraints

# Compile scenario with adjacency constraints
problem = compile_scenario(fm, scenario_name="adjacency_test")

# Add adjacency constraints
add_adjacency_constraints(problem, adj_matrix, dt_mapping)

# Set solver and solve
problem.solver("gurobi")
problem.solve()

# Check status and get solution
print(f"Solved: {problem.solved()}")
print(f"Objective value: {problem.z():.2f}")
solution = problem.solution()
for var_name, value in solution.items():
    if value > 0:
        print(f"  {var_name}: {value:.2f}")

Contiguous Area Requirements

Problem: Harvested areas must meet minimum size requirements for economic viability and regulatory compliance.

Mathematical Formulation:

For contiguous area requirements, we need to ensure that harvested areas form connected components of sufficient size:

\[\sum_{i \in C} x_{i,t} \geq A_{min} \cdot y_{C,t} \quad \forall C \in \mathcal{C}, \forall t\]

Where \(\mathcal{C}\) is the set of all possible connected components, \(A_{min}\) is the minimum contiguous area, and \(y_{C,t}\) is a binary variable indicating whether component \(C\) is harvested in period \(t\).

Implementation Approach:

  1. Identify Connected Components: Use graph theory to find connected components in the harvest schedule

  2. Size Constraints: Add constraints to ensure minimum size

  3. Binary Variables: Use binary variables to indicate which components are harvested

from scipy.sparse.csgraph import connected_components

def check_contiguous_areas(schedule, spatial_df, min_area):
    """Check if harvested areas meet contiguous area requirements.

    :param schedule: harvest schedule with area and period columns
    :param spatial_df: spatial inventory with geometry
    :param min_area: minimum contiguous area requirement
    :return: boolean indicating if requirements met
    """
    # Get harvested areas
    harvested = spatial_df[spatial_df['harvested'] == 1]

    if harvested.empty:
        return True

    # Find connected components
    n_components, labels = connected_components(
        adj_matrix[harvested.index.values],
        directed=False
    )

    # Check size of each component
    for comp_id in range(n_components):
        comp_area = harvested[labels == comp_id]['area_ha'].sum()
        if comp_area < min_area:
            print(f"Component {comp_id} too small: {comp_area:.2f} ha < {min_area} ha")
            return False

    return True

Raster-Based Spatial Analysis

Advantages of Raster Data:

  • Uniform resolution: Consistent cell sizes

  • Easy computation: Fast neighborhood operations

  • Integration: Compatible with GIS and remote sensing

  • Flexibility: Easy to add layers (elevation, slope, etc.)

Converting Vector to Raster:

import rasterio
from rasterio.features import geometry_mask

def vector_to_raster(geometry_df, raster_template, field='dt_code'):
    """Convert vector geometries to raster.

    :param geometry_df: geopandas DataFrame with geometries
    :param raster_template: rasterio Dataset for reference
    :param field: attribute field to use as values
    :return: raster array
    """
    # Create output array
    out_image = np.zeros(raster_template.shape, dtype=raster_template.dtypes[0])

    # Create mask from geometries
    shapes = [
        (geom, value)
        for geom, value in zip(geometry_df.geometry, geometry_df[field])
    ]

    # Rasterize
    out_image = rasterio.features.rasterize(
        shapes,
        out_shape=raster_template.shape,
        fill=0,
        transform=raster_template.transform,
        dtype=raster_template.dtypes[0]
    )

    return out_image

# Convert spatial data to raster
dt_raster = vector_to_raster(spatial_df, raster_template)
print(f"Raster shape: {dt_raster.shape}")

Using Raster in Optimization:

def add_raster_constraints(problem, raster, constraint_type='buffer'):
    """Add constraints based on raster data.

    :param problem: ws3 optimization problem
    :param raster: raster array
    :param constraint_type: type of constraint ('buffer', 'slope', etc.)
    """
    if constraint_type == 'buffer':
        # Get buffer zones (e.g., near water bodies)
        buffer_mask = raster == BUFFER_VALUE

        # Add constraints to prevent harvesting in buffer zones
        for row in range(raster.shape[0]):
            for col in range(raster.shape[1]):
                if buffer_mask[row, col]:
                    # Add constraint for this cell
                    cell_var = f"cell_{row}_{col}"
                    if cell_var in problem._vars:
                        problem.add_constraint(
                            name=f"buffer_{row}_{col}",
                            coeffs={cell_var: 1.0},
                            sense='=',
                            rhs=0.0
                        )

Case Study: Adjacency Constraints in TSA 24

Objective: Implement adjacency constraints for TSA 24 to ensure sustainable forest management.

Data Requirements:

  • Spatial inventory with polygon geometries

  • Adjacency relationships between polygons

  • Minimum harvest block size (e.g., 50 hectares)

Implementation Steps:

  1. Load spatial data and create adjacency matrix

  2. Compile base optimization scenario

  3. Add adjacency constraints

  4. Solve and compare results

  5. Analyze impact on harvest schedule

# Load TSA 24 spatial data
tsa24_spatial = gpd.read_file("data/tsa24_spatial.geojson")

# Create adjacency matrix
adj_matrix = create_adjacency_matrix(tsa24_spatial.geometry)

# Compile scenario
problem = compile_scenario(fm, scenario_name="tsa24_adjacency")

# Add adjacency constraints
add_adjacency_constraints(problem, adj_matrix, tsa24_spatial['dt_code'])

# Set solver and solve
problem.solver("gurobi")
problem.solve()

# Analyze results
solution = problem.solution()
for var_name, value in solution.items():
    if value > 0:
        print(f"  {var_name}: {value:.2f}")
print(f"Objective value: {problem.z():.2f}")

Summary

This chapter covered advanced spatial modeling techniques for forest optimization:

  • Adjacency constraints: Prevent harvesting adjacent areas

  • Contiguous area requirements: Ensure minimum harvest block sizes

  • Raster-based analysis: Use raster data for spatial constraints

  • Case study: Applied techniques to TSA 24

These techniques enable more realistic and practical forest management optimization that accounts for spatial structure and constraints.

Exercises

  1. Adjacency Matrix: Create an adjacency matrix for a simple 3x3 grid of forest polygons. Verify that each polygon is adjacent to its neighbors.

  2. Contiguous Areas: Modify the adjacency-constrained optimization to add contiguous area requirements. Test with different minimum area thresholds.

  3. Raster Constraints: Convert a slope raster to constraints that prevent harvesting on steep slopes (>30%).

  4. Performance Impact: Measure the impact of adjacency constraints on solve time and objective value. At what point do constraints make the problem infeasible?

  5. Real Data: Apply these techniques to a real forest inventory dataset. Compare results with and without spatial constraints.