Detailed code and steps used in the analysis

1 Environment Setup

Installed rasterio and gdal via pip within the Jupyter Notebook environment. Imported all required libraries for raster processing and visualization.

import numpy as np
import rasterio
import rasterio.plot
import matplotlib.pyplot as plt
import os

2 Reading Raster Data and Applying Suitability Criteria

Each input raster was opened using rasterio and evaluated against suitability thresholds to produce binary masks. Criteria included elevation between 300–700 m, road proximity, canopy height, area quality, hydrolines, cutblocks, MKRF bounds, and MGVI.

# Read DEM and apply elevation criterion
with rasterio.open("dem.tif") as dem_src:
    dem = dem_src.read(1).astype(float)
    profile = dem_src.profile
    dem[dem == dem_src.nodata] = np.nan

elevation_mask = (dem >= 300) & (dem <= 700)

# Read roads buffer and apply proximity criterion
with rasterio.open("RoadsBuff.tif") as roads_src:
    roads = roads_src.read(1)

roads_mask = (roads == 1)

# Combine all binary masks using logical AND
suitability = (elevation_mask & roads_mask
               & canopy_mask & area_mask
               & hydro_mask & cutblock_mask
               & mkrf_mask & mgvi_mask).astype("uint8")

3 Visualization — Multi-Panel Input Plot

A multi-panel matplotlib figure was created to display each input raster after applying suitability criteria. Yellow = suitable, purple = not suitable.

fig, axes = plt.subplots(2, 4, figsize=(16, 8))
fig.suptitle("Constraints for Maple Location", fontsize=14)

layers = {
    "Elevation": elevation_mask,
    "Roads": roads_mask,
    "Cutblocks": cutblock_mask,
    "Hydrolines": hydro_mask,
    "Area Quality": area_mask,
    "Canopy Height": canopy_mask,
    "MKRF Bounds": mkrf_mask,
    "MGVI": mgvi_mask
}

for ax, (title, layer) in zip(axes.flatten(), layers.items()):
    im = ax.imshow(layer, cmap="plasma")
    ax.set_title(title)

fig.colorbar(im, ax=axes[:, -1], label="Suitable")
plt.tight_layout()
plt.savefig("maple_constraints.png", dpi=150)
plt.show()

4 Visualization — Final Suitability Map

The final combined binary suitability raster was plotted using plt.imshow with a colorbar indicating suitable vs. not suitable areas.

plt.figure(figsize=(6, 10))
plt.imshow(suitability, cmap="plasma")
cbar = plt.colorbar(label="Suitable")
cbar.set_ticks([0, 1])
cbar.set_ticklabels(["0", "1"])
plt.title("Suitable Maple Plantation Locations")
plt.savefig("maple_suitability.png", dpi=150)
plt.show()

5 Exporting the Output Raster

The boolean output array was cast to uint8 format. The raster profile was updated and the output written to a GeoTiff file.

output_profile = profile.copy()
output_profile.update(dtype="uint8", count=1, nodata=0)

with rasterio.open("maple-farm-locations.tif", "w", **output_profile) as dst:
    dst.write(suitability, 1)

6 Overlay on True-Colour Satellite Image

The SPOT true-colour image was read and displayed using rasterio.plot.show. The suitability layer was masked so that non-suitable cells were transparent, then overlaid with 50% opacity.

with rasterio.open("Spot_data.tif") as spot_src:
    spot = spot_src.read()
    spot_transform = spot_src.transform

# Mask non-suitable cells
masked_suit = np.ma.masked_where(suitability == 0, suitability)

fig, ax = plt.subplots(figsize=(8, 14))
rasterio.plot.show(spot, transform=spot_transform, ax=ax)
ax.imshow(masked_suit, cmap="YlGn", alpha=0.5,
          extent=[...])  # extent from raster bounds
ax.set_title("Suitable Maple Plantation Locations")
plt.legend(handles=[plt.Rectangle((0,0),1,1, color="#9acd32", label="Viable Locations")])
plt.savefig("maple_overlay.png", dpi=150)
plt.show()