Detailed code and steps used in the analysis

1 Environment Setup

Imported required libraries and checked out the Spatial Analyst extension. Created a scratch geodatabase for intermediate outputs.

import arcpy
import pandas as pd
import os
import seaborn as sns
import matplotlib.pyplot as plt

# Check out Spatial Analyst extension
arcpy.CheckOutExtension("Spatial")

# Set up scratch geodatabase
scratch_gdb = r"C:\path\to\scratch.gdb"
if not arcpy.Exists(scratch_gdb):
    arcpy.management.CreateFileGDB(os.path.dirname(scratch_gdb),
                                   os.path.basename(scratch_gdb))
arcpy.env.workspace = scratch_gdb

2 Raster to Vector Conversion

Converted the Murray_Mangrove_Loss_Year.tif raster to a polygon feature class using arcpy.conversion.RasterToPolygon. Used NO_SIMPLIFY to preserve geometry and mapped the Value field to retain loss year values.

input_raster = r"C:\path\to\Murray_Mangrove_Loss_Year.tif"
output_polygons = os.path.join(scratch_gdb, "mangrove_loss_polygons")

arcpy.conversion.RasterToPolygon(
    in_raster=input_raster,
    out_polygon_features=output_polygons,
    simplify="NO_SIMPLIFY",
    raster_field="Value"
)

3 Data Cleaning

Renamed the gridcode field to loss_year and removed all zero-value (no-data) polygons using an UpdateCursor.

# Rename gridcode to loss_year
fields = [f.name for f in arcpy.ListFields(output_polygons)]
if "gridcode" in fields:
    arcpy.management.AlterField(output_polygons, "gridcode",
                                 "loss_year", "loss_year")

# Remove zero-value records
with arcpy.da.UpdateCursor(output_polygons, ["loss_year"]) as cursor:
    for row in cursor:
        if row[0] == 0:
            cursor.deleteRow()

# Validate removal
count = sum(1 for row in arcpy.da.SearchCursor(
    output_polygons, ["loss_year"]) if row[0] == 0)
print(f"Zero-value records remaining: {count}")

4 Temporal Trend Visualization

Extracted the attribute table into a pandas DataFrame and generated a seaborn histogram to visualize deforestation trends over time.

# Extract to DataFrame
df = pd.DataFrame.spatial.from_featureclass(output_polygons)

# Plot histogram
plt.figure(figsize=(10, 5))
sns.histplot(data=df, x="loss_year", bins=20, color="#2e7d32")
plt.title("Mangrove Loss Year — Ambanja and Ambaro Bays")
plt.xlabel("Year of Loss")
plt.ylabel("Count")
plt.tight_layout()
plt.savefig("mangrove_histogram.png", dpi=150)
plt.show()

5 Quantifying Loss by Administrative Unit

Copied commune and fokontany boundaries from the AAB geodatabase, performed spatial intersections, calculated area in km², and summarized by administrative unit.

aab_gdb = r"C:\path\to\AAB.gdb"

for admin_layer, name_field in [("commune", "COMMUNE"), ("fokontany", "FOKONTANY")]:
    src = os.path.join(aab_gdb, admin_layer)
    dst = os.path.join(scratch_gdb, admin_layer)
    arcpy.management.CopyFeatures(src, dst)

    # Spatial intersection
    intersect_out = os.path.join(scratch_gdb, f"loss_{admin_layer}")
    arcpy.analysis.Intersect([output_polygons, dst], intersect_out)

    # Calculate area in km²
    arcpy.management.CalculateGeometryAttributes(
        intersect_out, [["area_km2", "AREA"]], area_unit="SQUARE_KILOMETERS")

    # Summarize
    stats_out = os.path.join(scratch_gdb, f"stats_{admin_layer}")
    arcpy.analysis.Statistics(
        intersect_out, stats_out,
        [["area_km2", "SUM"], ["loss_year", "MEAN"]],
        case_field=name_field)

    # Export to Excel
    arcpy.conversion.TableToExcel(stats_out, f"loss_by_{admin_layer}.xlsx")

6 Buffer Analysis — Proximity to Paths and Rivers

Created multiple ring buffers at 1,000 m, 2,000 m, and 5,000 m around permanent rivers and paths, then intersected with mangrove loss polygons to quantify loss per distance band.

for feature, label in [("permanent_rivers", "rivers"), ("paths", "paths")]:
    src = os.path.join(aab_gdb, feature)
    dst = os.path.join(scratch_gdb, feature)
    arcpy.management.CopyFeatures(src, dst)

    # Multiple ring buffers
    buffer_out = os.path.join(scratch_gdb, f"buffer_{label}")
    arcpy.analysis.MultipleRingBuffer(
        dst, buffer_out,
        distances=[1000, 2000, 5000],
        distance_unit="Meters",
        buffer_unit="Meters",
        outside_polygons_only="FULL")

    # Intersect with loss polygons
    loss_buffer = os.path.join(scratch_gdb, f"loss_buffer_{label}")
    arcpy.analysis.Intersect([output_polygons, buffer_out], loss_buffer)

    # Calculate area
    arcpy.management.CalculateGeometryAttributes(
        loss_buffer, [["area_km2", "AREA"]], area_unit="SQUARE_KILOMETERS")

    # Summarize by buffer ring
    stats_out = os.path.join(scratch_gdb, f"stats_buffer_{label}")
    arcpy.analysis.Statistics(
        loss_buffer, stats_out,
        [["area_km2", "SUM"]], case_field="distance")