Protected Area Data Explained

On this page

Problem statement

"How much of this is protected?" sounds like a spatial join and is not. Protected area datasets contain overlapping designations, marine and terrestrial parts of the same site, points standing in for polygons that were never digitised, and categories that mean entirely different things — a strict nature reserve and a landscape where farming continues are both in the same layer.

British Columbia publishes 929 polygons under four designations covering 111,299 km². Whether that is 12% of the province depends on questions the layer does not answer for you.

Quick answer

Load the polygons, look at what the designations are, and dissolve before measuring:

import geopandas as gpd
from pa import protected_areas          # the service wrapper below

pas = protected_areas()                 # BC Albers, EPSG:3005
pas["km2"] = pas.area / 1e6

print(f"{len(pas)} polygons, {pas.km2.sum():,.0f} km² summed")
print(pas.PROTECTED_LANDS_DESIGNATION.value_counts().to_string())

dissolved = pas.union_all()
print(f"\nsummed area:   {pas.km2.sum():>10,.0f} km²")
print(f"dissolved area:{dissolved.area / 1e6:>10,.0f} km²")
print(f"difference:    {pas.km2.sum() - dissolved.area / 1e6:>10,.0f} km² of overlap")
929 polygons, 111,299 km² summed
PROTECTED_LANDS_DESIGNATION
PROVINCIAL PARK       694
ECOLOGICAL RESERVE    148
PROTECTED AREA         85
RECREATION AREA         2

summed area:      111,299 km²
dissolved area:   111,217 km²
difference:            81 km² of overlap

The summed figure is what a spreadsheet gives you; the dissolved figure is the area of ground that is protected. Here they differ by only 81 km² — this provincial layer is almost non-overlapping, which is worth knowing rather than assuming. In the global database, where a national park, a biosphere reserve and a World Heritage site can cover the same valley, the gap is percentage points rather than fractions of one.

Protected area designations in one jurisdiction with counts and areas. four designations, one layerdesignationsiteswhat it usually permitsprovincial park694recreation, some infrastructureecological reserve148research access onlyprotected area85varies by order in councilrecreation area2the most permissive of the fourecological reserves are 16% of the sites and 1.5% of the protected areaBritish Columbia: 929 polygons, 111,217 km² dissolved spatialworkflow.io
Four designations in one layer; they do not mean the same thing on the ground.

What the layer contains, designations, dissolving, equal-area

1. Find out what the layer actually contains

Terrestrial, marine, proposed, de-gazetted, private conservancies: the mix varies by publisher and is rarely in the layer's name.

2. Read the designation, not just the boundary

IUCN categories I to VI run from strict reserve to sustainable use. A map that treats them alike answers a question nobody asked.

3. Dissolve before measuring area

Overlapping designations are normal — a park inside a biosphere reserve inside a World Heritage site. Summing their areas counts the same ground repeatedly.

4. Measure in an equal-area projection

Protected area statistics are area statistics; a geographic CRS gives square degrees, and Web Mercator inflates high latitudes enormously.

5. Check the denominator

"12% protected" needs the total area it is a percentage of — land only, land and sea, or the administrative area including water.

6. Note what is missing

Sites recorded only as a point with a reported area, sites whose boundaries are restricted, and recent designations not yet in the layer.

Overlapping protected area designations and the difference between summed and dissolved area. summing designations counts shared ground twicenational parkbiosphere reservecounted twicesummed111,299 km²dissolved111,217 km²in this provincial layer the overlap is only 81 km²; in the global databaseit is percentage points, because designations routinely nest spatialworkflow.io
Summing designations counts shared ground once per designation.

British Columbia's protected areas measured

The data used on this page

The protected areas are British Columbia's own layer, BC Parks, Ecological Reserves and Protected Areas, from the province's ArcGIS feature service in BC Albers (EPSG:3005). pa.py wraps the download so every example starts from the same cached file.

import os

with open("pa.py", "w") as f:
    f.write('''
import os
import geopandas as gpd
import requests

URL = ("https://services6.arcgis.com/ubm4tcTYICKBpist/arcgis/rest/services/"
       "British_Columbia_Parks_Ecological_Reserves_and_Protected_Areas/FeatureServer/0/query"
       "?where=1%3D1&outFields=PROTECTED_LANDS_NAME,PROTECTED_LANDS_DESIGNATION&outSR=3005&f=geojson")


def protected_areas():
    if not os.path.exists("bc_protected.geojson"):             # cache the service's own GeoJSON
        with open("bc_protected.geojson", "wb") as f:
            f.write(requests.get(URL, timeout=300).content)
    return gpd.read_file("bc_protected.geojson")
''')

from pa import protected_areas

print(len(protected_areas()), "polygons")
929 polygons

Example 1 — the size distribution, which decides everything

import numpy as np

print(f"median size {pas.km2.median():.2f} km², mean {pas.km2.mean():.1f} km²")
for lo, hi, label in [(0, 1, "under 1 km²"), (1, 10, "1–10 km²"),
                      (10, 100, "10–100 km²"), (100, 1_000, "100–1,000 km²"),
                      (1_000, 1e9, "over 1,000 km²")]:
    sel = (pas.km2 >= lo) & (pas.km2 < hi)
    print(f"  {label:<16} {sel.sum():>4} sites, {pas.km2[sel].sum():>9,.0f} km² "
          f"({pas.km2[sel].sum() / pas.km2.sum():>5.1%} of the area)")
print(f"\nthe ten largest sites hold {pas.km2.nlargest(10).sum() / pas.km2.sum():.0%} "
      f"of the protected area")
median size 2.62 km², mean 119.8 km²
  under 1 km²       310 sites,       113 km² ( 0.1% of the area)
  1–10 km²          327 sites,     1,187 km² ( 1.1% of the area)
  10–100 km²        170 sites,     5,950 km² ( 5.3% of the area)
  100–1,000 km²      98 sites,    31,312 km² (28.1% of the area)
  over 1,000 km²     24 sites,    72,736 km² (65.4% of the area)

the ten largest sites hold 47% of the protected area

A protected area network is nearly always a few enormous sites and a long tail of small ones. Any statistic that weights sites equally — mean size, percentage of sites meeting a criterion — is dominated by the tail, and any statistic that weights by area is dominated by a handful of places.

Example 2 — what a percentage needs

land_km2 = 944_735            # British Columbia's land area, published figure
total_km2 = 944_735 + 39_000  # plus fresh water, roughly

protected = dissolved.area / 1e6
print(f"protected area (dissolved): {protected:,.0f} km²")
print(f"  as a share of land:            {protected / land_km2:.1%}")
print(f"  as a share of land plus water: {protected / total_km2:.1%}")

strict = pas[pas.PROTECTED_LANDS_DESIGNATION == "ECOLOGICAL RESERVE"]
print(f"\nstrict reserves only: {len(strict)} sites, "
      f"{strict.area.sum() / 1e6:,.0f} km² ({strict.area.sum() / 1e6 / land_km2:.2%} of land)")
print(f"  median strict reserve size: {strict.area.median() / 1e6:.2f} km²")
protected area (dissolved): 111,217 km²
  as a share of land:            11.8%
  as a share of land plus water: 11.3%

strict reserves only: 148 sites, 1,628 km² (0.17% of land)
  median strict reserve size: 1.61 km²

The same layer supports "protected" figures an order of magnitude apart, depending on which designations count and what the denominator is. Both are defensible; publishing one without saying which is not.

Example 3 — validity, before any overlay

from shapely.validation import make_valid

invalid = ~pas.geometry.is_valid
print(f"invalid geometries: {invalid.sum()} of {len(pas)}")
if invalid.any():
    print(pas.loc[invalid, ["PROTECTED_LANDS_NAME", "km2"]].head(3).to_string(index=False))

fixed = pas.copy()
fixed["geometry"] = [make_valid(g) for g in fixed.geometry]
print(f"after make_valid: {(~fixed.geometry.is_valid).sum()} invalid")
print(f"area change from repair: "
      f"{abs(fixed.area.sum() - pas.area.sum()) / 1e6:,.3f} km²")

multi = pas.geometry.geom_type.value_counts()
print(f"\ngeometry types: {multi.to_dict()}")
print(f"parts in the largest multipart site: "
      f"{max(len(g.geoms) if g.geom_type == 'MultiPolygon' else 1 for g in pas.geometry)}")
invalid geometries: 40 of 929
          PROTECTED_LANDS_NAME         km2
               ARROWSTONE PARK   61.715546
ATLIN/A TEIX'GI AAN TLEIN PARK 2299.288964
             BROMLEY ROCK PARK    1.508038
after make_valid: 0 invalid
area change from repair: 95.678 km²

geometry types: {'Polygon': 619, 'MultiPolygon': 310}
parts in the largest multipart site: 265

Forty of the 929 polygons are invalid, and repairing them changes the total area by 96 km² — more than the 81 km² of overlap the dissolve removed. Both numbers are small against 111,000 km², and both are larger than the difference between several published "percentage protected" figures for this province. The repair has to happen before the overlay either way, because an invalid polygon makes union_all and every spatial predicate unreliable.

Why overlapping designations are the norm

Why overlapping designations are the norm

Protection accretes. A park is designated, part of it becomes a wilderness zone, the whole of it is later inside a UNESCO biosphere reserve, and a marine portion is added under separate legislation. Each layer is correct; the union is what is protected, and the sum is a number with no referent.

Why the category matters more than the boundary

IUCN category I permits almost no human use; category VI explicitly allows sustainable resource extraction. A "protected area" map that merges them supports conclusions about neither, and the gap is not small — in many countries the majority of the protected estate is in the permissive categories.

Why the point records exist

The world database records sites whose boundaries were never digitised or cannot be released as a point with a reported area. They are legitimate entries and cannot be used in an overlay; dropping them silently understates coverage, and buffering them to their reported area invents a shape.

Protected area count and area by size class. a few huge sites and a long tailover 1,000 km² (24 sites)65.4% of the area100–1,000 km² (98)28.1%10–100 km² (170)5.3%1–10 km² (327)1.1%under 1 km² (310)0.1% — a third of the sites037.5k km²75k km²Median site 2.6 km², mean 119.8 km²: any per-site statistic is dominated by the tail spatialworkflow.io
24 of 929 sites hold two thirds of the area; the other 905 are the tail.

WDPA, marine sites and OECMs

  • WDPA is the global database; its use requires citing the monthly release, because it changes every month.
  • Marine sites are in the same layer and must be excluded from terrestrial statistics.
  • Proposed and de-gazetted sites carry a status field; including them changes the answer.
  • OECMs — other effective area-based conservation measures — are a separate, growing layer.
  • Indigenous protected areas are recorded inconsistently and are often missing entirely.
  • Reported area and GIS area disagree, exactly as with any other polygon dataset.

FAQ

What is a protected area dataset?

A layer of designated conservation sites with their boundaries, designations and status — WDPA globally, plus national and provincial layers.

Why can't I just sum the areas?

Because designations overlap. Dissolve the geometry first, or you count the same ground more than once.

What do the IUCN categories mean?

They range from strict nature reserve (I) to sustainable use of natural resources (VI); they are not interchangeable.

Why are some sites points?

Their boundaries were never digitised or cannot be released. They carry a reported area and cannot be used in an overlay.

Which denominator should a percentage use?

Say which: land area, land and inland water, or total administrative area. They differ by several per cent.

Does the data change?

Constantly. WDPA publishes monthly; national layers update as sites are designated. Cite the release you used.