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.
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.
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.
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.
Internal links
- How to find which species occur in protected areas — the overlay this prepares for
- Biodiversity metrics explained — what gets summarised per site
- Species occurrence data explained — the records joined to these polygons
- How to dissolve polygons in Python (GeoPandas) — the operation that fixes the double count
- How to calculate area and distance in GeoPandas (correctly) — the projection question
- Exposure counts are double-counted or missed — the same arithmetic in hazard work
- How to compare forest loss inside and outside protected areas — signed distance from the boundary
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.