Spatial Analysis for
Wildlife Conservation

The Big Picture

Imagine you have years of wildlife survey data — GPS points from camera traps, ranger observations, interview records, and digitized maps. Each piece tells a small story. Our pipeline stitches those stories together into a single, spatially explicit picture of biodiversity across the AOI (Area of Interest): roughly 110 km² around the Upper Cisokan Pumped Storage Hydroelectric Power Plant (REEPS) in West Java.

The Pipeline

Meet the Focal Species

Eleven species — mostly IUCN-listed — form the core of our analysis. This is the cast of characters you will track through every module.

Synthetic Data Generator

Before touching real data, let's generate a synthetic dataset with realistic parameters. Click the button to produce ~400 random occurrence records.

Showing first 15 of records. Full dataset stored in memory.

Parsing KMZ Files

A KMZ file is simply a zipped archive containing one or more KML files. KML (Keyhole Markup Language) is an XML format that stores geographic features — points, lines, polygons — along with metadata. Our surveys often export sighting locations as KMZ from Google Earth.

Code ↔ English

Python
with zipfile.ZipFile(kmz_path) as z:
    for name in z.namelist():
        if name.endswith('.kml'):
            content = z.open(name).read().decode('utf-8')
            root = ET.fromstring(content)
            for pm in root.iter('{http://www.opengis.net/kml/2.2}Placemark'):
                point = pm.find('.//{http://www.opengis.net/kml/2.2}coordinates')
                coords = point.text.strip().split(',')
                lon, lat = float(coords[0]), float(coords[1])

Species Harmonization

The same animal may appear under different names across datasets: local Indonesian names, outdated scientific names, or variant spellings. Species harmonization maps all variants to a single canonical scientific name.

Python
SPECIES_MAP = {
    'Kukang': 'Nycticebus javanicus',
    'Lutung': 'Trachypithecus auratus',
    'Owa Jawa': 'Hylobates moloch',
    'Kucing Hutan': 'Prionailurus bengalensis',
    'Trenggiling': 'Manis javanica',
}
df['Species'] = df['Species'].map(lambda x: SPECIES_MAP.get(x, x))

Interactive: Match Local Names to Scientific Names

Drag each local Indonesian name to its correct scientific name.

Why Hexagons?

Before we grid our study area, we need to pick a shape. Squares seem obvious, but hexagons are the gold standard for ecological spatial analysis. Here is why:

H3 in Action

Uber's H3 library assigns every point on Earth to a hexagonal cell at a chosen resolution. Here is how we apply it to our wildlife records:

Code ↔ English

Python
import h3

# Assign each wildlife sighting to its
# enclosing hexagonal cell
df['h3_index'] = df.apply(
    lambda r: h3.latlng_to_cell(
        r['Latitude'],
        r['Longitude'], 8), axis=1)

# Convert H3 cell to polygon for mapping
boundary = h3.cell_to_boundary(cell)
polygon = Polygon(
    [(lng, lat) for lat, lng in boundary])

Resolution Matters

H3 supports 16 resolution levels (0–15). Too coarse and you mix habitats; too fine and most cells are empty. The sweet spot depends on your data density and the organisms' home ranges.

Interactive: Explore Hex Neighbors

Click any hexagon to select it. Its K1 ring (immediate neighbors, like h3.grid_disk(cell, 1)) lights up in green, and the K2 ring (neighbors-of-neighbors) appears in lighter green. This is the basis for spatial autocorrelation analysis.

The Diversity Question

Consider two hexagonal cells in our study area. Which one is more biodiverse?

Shannon Entropy (H′)

Borrowed from information theory, Shannon entropy measures the "surprise" or uncertainty in a community. Higher H′ means more species and/or more even distribution — a key measure of alpha diversity.

Code ↔ English

Python
from math import log
from collections import Counter

counts = Counter(cell_species)
# e.g., {'Leopard': 10, 'Gibbon': 5,
#        'Pangolin': 3}
N = sum(counts.values())  # total = 18

H = -sum(
    (n/N) * log(n/N)
    for n in counts.values()
    if n > 0
)

Simpson's Index (D) & Pielou Evenness (J)

Two complementary metrics that capture different facets of diversity. Simpson asks "what is the probability two random individuals are different species?" while Pielou normalizes Shannon by its theoretical maximum.

# Probability that two random
# individuals are different species
D = 1 - sum(
    (n/N) ** 2
    for n in counts.values()
)

# Shannon normalized by max
# possible Shannon (ln S)
J = H / log(S) if S > 1 else 0

Interactive: Diversity Calculator

Adjust the sliders to change species counts and watch H′, D, and J update in real time. Try making one species dominant, then redistribute evenly.

Chao1 — How Many Species Are Hiding?

Observed species richness is always an undercount. Chao1 uses the frequency of rare species — singletons (species seen once) and doubletons (seen twice) — to estimate the true species pool.

Code ↔ English

Python
S_obs = len(counts)  # observed richness

# Count singletons and doubletons
f1 = sum(1 for v in counts.values()
         if v == 1)  # singletons
f2 = sum(1 for v in counts.values()
         if v == 2)  # doubletons

if f2 > 0:
    chao1 = S_obs + (f1 ** 2) / (2 * f2)
elif f1 > 0:
    chao1 = S_obs + f1 * (f1 - 1) / 2

completeness = S_obs / chao1 * 100

Why Not Just Count Species?

Imagine two cells. Cell A has 3 common species — Leopard Cat, Palm Civet, Mouse-deer — all Least Concern. Cell B has 2 critically endangered species — Javan Leopard and Sunda Pangolin. Which cell matters more for conservation?

Five CPI Components

PCA (Principal Component Analysis) revealed which dimensions carry the most variance. The resulting weights tell us how much each component matters in the final score.

Min-Max Normalization

Before combining scores, each component must be on the same 0–1 scale. Otherwise, a metric measured in hundreds would drown out one measured in single digits.

Code ↔ English

Python
for col in ['Richness', 'Diversity',
            'Connectivity', 'Co_occurrence',
            'Threatened']:
    vmin, vmax = df[col].min(), df[col].max()
    if vmax > vmin:
        df[f'n_{col}'] = (df[col] - vmin) / (vmax - vmin)
    else:
        df[f'n_{col}'] = 0

Weighted Sum

Once each component is normalized, we multiply by its PCA-derived weight and add them up. The result: one number per cell.

Code ↔ English

Python
weights = {
    'Richness': 0.28,
    'Diversity': 0.18,
    'Connectivity': 0.22,
    'Co_occurrence': 0.12,
    'Threatened': 0.20
}

df['Priority_Index'] = sum(
    df[f'n_{k}'] * v
    for k, v in weights.items()
)

Priority Tiers

The continuous CPI score is binned into four actionable tiers, each with a clear management implication.

Interactive: CPI Calculator

Adjust each component on a 0–10 scale. The calculator applies min-max normalization (treating 0 as min, 10 as max) and the PCA weights to produce a CPI score and tier classification.

The Airport Hub Metaphor

Think of occupied hex cells as airports. Each airport has flights (edges) to neighboring airports. Empty cells between clusters are like missing routes. Some empty cells — if "built" — could connect two isolated airport networks. That empty cell is a stepping stone.

Building the Network

We use NetworkX, Python's graph theory library, to turn the hex grid into a network where each cell is a node and touching cells share an edge.

Code ↔ English

Python
import networkx as nx

G = nx.Graph()
for cell in all_h3:
    G.add_node(cell)
for cell in all_h3:
    for nbr in h3.grid_disk(cell, 1):
        if nbr != cell and nbr in G:
            G.add_edge(cell, nbr)

Connectivity Metrics

Each cell's position in the network reveals its ecological role. Interior cells are safe havens; keystone cells are structurally critical.

Stepping Stone Analysis

Empty cells are scored by how effectively they could bridge gaps between occupied habitat. Betweenness centrality measures how many shortest paths pass through a node — high betweenness means the cell is a critical bridge.

Code ↔ English

Python
betweenness = nx.betweenness_centrality(G)

score = (
    len(k1_occ) * 3
    + k2_richness * 0.5
    + betweenness.get(cell, 0) * 10
)
tier = ('High' if score > 5
        else ('Medium' if score > 2
        else 'Low'))

Resistance Surface

Not all land is equally permeable. A resistance surface assigns a movement cost to each land cover type. Low resistance = easy passage; high resistance = major barrier.

Interactive: Stepping Stone Explorer

This simplified 7×5 hex grid shows occupied cells (green) and empty cells (gray). Click or tap any empty cell to calculate its stepping stone score based on occupied neighbors and estimated betweenness. The cell will be colored by its corridor tier.

Why Interactive Maps?

Static images are fine for a report, but they cannot zoom into a single hex cell, toggle layers on and off, or show a tooltip when you hover. Folium wraps the powerful Leaflet.js mapping library in pure Python — one function call produces a self-contained .html file you can open in any browser.

Building a Base Map

Every Folium map starts with a center point and a basemap tile layer.

Code ↔ English

Python
import folium
import branca.colormap as cm

center = [h3_all.geometry.centroid.y.mean(),
          h3_all.geometry.centroid.x.mean()]
m = folium.Map(location=center,
               zoom_start=13,
               tiles='CartoDB positron')

Adding Hex Layers with Styling

Each hex cell is drawn as a GeoJSON polygon, colored by its metric value. Hovering reveals a tooltip with the cell’s details.

Code ↔ English

Python
cmap_r = cm.LinearColormap(
    ['#ffffcc', '#fd8d3c', '#bd0026'],
    vmin=0, vmax=max_r,
    caption='Species Richness')

folium.GeoJson(
    h3_div.__geo_interface__,
    name='Species Richness',
    style_function=lambda f: {
        'fillColor': cmap_r(
            f['properties'].get('Richness__S', 0)),
        'color': '#333333',
        'weight': 1,
        'fillOpacity': 0.6},
    tooltip=folium.GeoJsonTooltip(
        fields=['h3_index', 'Richness__S',
                'Shannon__H'],
        aliases=['Cell', 'Richness',
                 "Shannon H'"])
).add_to(m)

Layer Control

Folium’s built-in LayerControl adds a toggle panel to the map. Each GeoJSON overlay gets a checkbox — users can switch between diversity, priority, and corridor layers without reloading.

The 7 Maps Produced

The pipeline generates seven interactive HTML maps, each focused on a different analytical theme.

Why Static Figures?

Interactive maps cannot be printed at 300 DPI. Peer reviewers expect precise axis labels, scale bars, and north arrows. Matplotlib produces publication-quality vector graphics (PDF, SVG) and high-resolution raster images (PNG) — exactly what journals require.

Loading a Satellite Basemap

Our figures use a Sentinel-2 satellite image as a basemap. Rasterio reads the GeoTIFF and reprojects it from UTM to WGS 84 (EPSG:4326).

Code ↔ English

Python
import rasterio
from rasterio.warp import (
    reproject, Resampling,
    calculate_default_transform)

with rasterio.open(S2_TIF) as src:
    transform, width, height = \
        calculate_default_transform(
            src.crs, "EPSG:4326",
            src.width, src.height,
            *src.bounds)
    rgb = np.zeros((3, height, width),
                   dtype=np.uint8)
    for band in range(1, 4):
        reproject(
            source=rasterio.band(src, band),
            destination=rgb[band - 1],
            dst_transform=transform,
            dst_crs="EPSG:4326",
            resampling=Resampling.bilinear)

Building the Figure

With the satellite image loaded, we layer hex polygons on top and save at print resolution.

Code ↔ English

Python
fig = plt.figure(figsize=(7.0, 5.4),
                 facecolor="white")
ax = fig.add_axes([0.02, 0.02,
                   0.78, 0.96])

ax.imshow(s2_rgba, extent=s2_extent,
         origin="upper",
         interpolation="bilinear")

gdf_unocc.plot(ax=ax, color="#E8E8E8",
    alpha=0.22, edgecolor="#555",
    linewidth=0.35)

gdf_occ.plot(ax=ax, color=colors,
    alpha=0.65, edgecolor="#555",
    linewidth=0.35)

fig.savefig("richness.png",
    dpi=300, bbox_inches="tight")
fig.savefig("richness.pdf",
    bbox_inches="tight")

Figure Gallery

The pipeline produces seven figure sets — each with multiple panels covering different metrics and composite layouts.

Publication Styling Tips

Small details separate a good figure from a publishable one. Follow these steps for every map figure.