Introduction
I’ve always had a deep appreciation for forests and the role they play in maintaining our planet’s health. But as reforestation efforts become more critical in the fight against climate change and biodiversity loss, I began asking myself a question: Where should we focus these efforts for the greatest impact? That question sparked the beginning of this project—a journey into the world of geospatial analysis using Python and ArcGIS.
This wasn’t just about generating a final suitability map. It was also a personal challenge to grow my programming skills as a beginner Python user, and to explore how scripting can enhance the power of traditional GIS workflows. I hope sharing my experience can help others who are also trying to bridge that gap between GIS and Python.
Project Goal
My main objective was to develop a reforestation suitability model: a raster map highlighting areas within my study region that would be most suitable for reforestation. I focused on the HU4 watershed district spanning parts of Vermont and New Hampshire. I selected this part of the country because it was a different location to past projects.
To answer the question of where reforestation efforts would likely have the greatest impact, I used multiple environmental datasets and Python libraries to process, analyze, and combine spatial layers based on criteria such as slope, precipitation, land cover, temperature, and proximity to water.
Gathering & Preparing the Data
I collected datasets from a variety of sources:
- Land Cover (NLCD) from USGS
- Slope from a merged DEM (Digital Elevation Model)
- Precipitation and Temperature from WorldClim
- Water Proximity using OpenStreetMap and HydroSHEDS
Each dataset had its own quirks—different coordinate systems, spatial resolutions, and file formats. Some were multi-part shapefiles, others were monthly raster layers or global datasets that needed clipping. Learning how to clip, reproject, and align everything using Python was a crash course in geospatial data handling.
Highlights included:
- Jumping through various hoops to obtain datasets, such as submitting a data request for the area of study
- Merging DEM tiles and generating a slope raster using Rasterio
- Rasterizing OSM water features and calculating Euclidean distance to water
- Reprojecting and clipping all layers to match my HU4 watershed boundary
Defining Suitability Criteria
I created a set of rules to guide reforestation site selection. Here’s what I looked for:
| Layer | Suitability Logic |
|---|---|
| Slope | Flat to moderate slopes preferred (<20 degrees) |
| Land Cover | Avoid already-existing forests and developed areas; prefer shrubland |
| Precipitation | Higher rainfall (>1000 mm/year) = more suitable |
| Temperature | Moderate climates (~5–15°C) ideal |
| Water Proximity | Closer to water bodies = better |
Each layer was reclassified to a 0–1 scale using Python. For example, slopes under 20 degrees were assigned a 1.0 (ideal), while steep slopes over 30 got a 0.1. Land cover types like shrubland and grassland were considered suitable, while existing forests and urban areas were masked out.
Reclassification & Weighting
Once reclassified, I normalized all rasters to the same spatial resolution and extent using Rasterio’s reproject() function. I then applied weights to reflect the relative importance of each factor:
- Slope: 0.25
- Land Cover: 0.30
- Precipitation: 0.15
- Temperature: 0.15
- Water Proximity: 0.15
Combining them was straightforward math: multiply each raster by its weight, then sum the results to get a final suitability score for each pixel.
The Final Suitability Map
The result was a raster called Reforestation_Suitability.tif, with values ranging from 0 (least suitable) to 1 (most suitable). It was incredibly satisfying to open that final file in ArcGIS and see the areas where all the conditions best aligned, represented by darker shades of green on the final map at the top of this post.
This map can now be used as the basis for further analysis—like applying a threshold to isolate only the top 10% of suitable areas, or overlaying with protected lands or ownership parcels.
Reflections
This project taught me far more than just the basic Python syntax one might learn about in an online course. It taught me how to troubleshoot unexpected errors, how to think through spatial logic, and how to persevere when my script didn’t work the first (or fifth) time. I now feel much more confident working with Rasterio, GeoPandas, and large-scale environmental datasets.
I also realized that scripting opens up a whole new level of efficiency and repeatability in GIS. What would’ve taken hours of manual clicking & typing can now be automated with just a few lines of Python code.
And that is why Python is such a valuable skill to have in the field of GIS.