ENCN205 — Demo Answer Guide

• Coastline — LineString (line tracing the shoreline)
• Roads — Line (centrelines of road corridors)
• Elevation Contours — LineString (lines at elevation intervals, not filled areas)
• Buildings — Polygon (footprint outlines of structures)
• Marae — Point (location markers)
• Structures — Point (Transpower power pylons)

You need at minimum two layers:

• Flood zone / Elevation contours — to define the area that would be flooded at 2m sea level rise
• Buildings — to identify which structures fall within (intersect) the flood zone

Optionally, you might also want the coastline to orient yourself and understand where the sea is relative to the land. This is a fundamental GIS operation: overlaying two layers to find spatial relationships.

The flood extent grows significantly between 1m and 3m, but the relationship is not linear. At 1m, only a very narrow coastal strip is affected (~0.06 km²). At 2m, the area jumps dramatically (~2.4 km²) because the low-lying coastal flat behind the dunes becomes inundated. At 3m, it expands further (~3.3 km²) but the rate of increase slows because the land rises more steeply inland.

This is because coastal topography is not uniform — there's a flat zone near the coast (the estuary/dune system) and rising ground further inland. The relationship between sea level and flooded area depends on the shape of the terrain (hypsometric curve).

Any three of:

• Flood zone / Elevation — defines which areas will be inundated and which are safe (determines who needs evacuating and where they can go)
• Roads — the actual routes people would use to evacuate; need to check which roads remain passable above flood level
• Buildings — identifies where people are (residential density) and therefore where evacuation demand is highest
• Marae — potential evacuation centres / welfare facilities; culturally significant sites that may need protection
• Structures — infrastructure like bridges and sea walls that affect route viability

Missing layers could include:

• Population / census data — meshblock population counts to estimate how many people are affected
• High-resolution DEM (raster) — the contours are coarse; an 8m or 1m LiDAR DEM would give much more accurate flood extents
• Stormwater / drainage network — pipes and drains that affect where water goes during flooding
• Emergency services — hospitals, fire stations, police stations; proximity to affected areas
• Land use zoning — residential vs commercial vs industrial affects vulnerability and response priorities
• Historical flood extent data — records of past flood events for calibration
• Property values — for economic loss estimation

When elevation is rounded to the nearest 5m, the flood zone boundary becomes much coarser and blocky — it either massively overestimates flooding (capturing large areas rounded down to 0m) or underestimates it (areas rounded up to 5m escape the threshold). The boundary loses its detailed shape and becomes stepped.

This demonstrates that input data quality directly controls output accuracy. A 5m-precision DEM cannot distinguish between 1m and 4m elevation — critical for flood modelling where 1-2m differences determine which streets flood.

At around 50-100m positional error, buildings near the flood zone boundary start being misclassified — buildings actually in the flood zone may appear safe, and vice versa. For resource consent, this absolutely matters: a building consent might be wrongly granted for a flood-prone site, or wrongly denied for a safe site.

In NZ, resource consents under the RMA require reliable spatial data. Consumer-grade GPS (~5-10m accuracy) is usually sufficient for regional planning, but site-specific decisions need survey-grade accuracy (~0.01m).

Marae are culturally significant sites for Māori communities. Publishing precise locations could:

• Enable unwanted visits or cultural tourism without appropriate tikanga (protocols)
• Expose sites to vandalism or theft of taonga (treasures)
• Remove iwi/hapū control over who accesses information about their spaces

Te Mana Raraunga asserts that Māori have inherent rights over data about Māori people, places, and taonga. This means data about marae should be managed with the consent and governance of the relevant iwi/hapū, not simply published openly because it's "public data." The NZ Privacy Act 2020 provides a legal framework, but Māori data sovereignty goes further — it's about rangatiratanga (self-determination) over data.

This dataset is likely NOT fit for purpose. Three concerns:

• Temporal currency — 20-year-old data doesn't reflect the 2010/2011 Canterbury earthquakes which changed ground levels by up to 0.5m in some coastal areas. The land surface has literally moved.
• Sea level change — Since 2005, sea level has risen ~3-4cm and projections have been updated. The flood extents are based on outdated climate scenarios.
• Land use change — 20 years of development has changed impervious surfaces, drainage patterns, and building stock. The analysis context has changed even if elevation hasn't.

Additional: resolution (1m is good for this purpose) and datum (NZGD2000 is correct) are acceptable, but the temporal issues are deal-breakers.

Using the root-sum-of-squares error propagation formula:

Combined uncertainty = √(0.5² + 0.3²) = √(0.25 + 0.09) = √0.34 ≈ ±0.58m

This means the flood boundary could be off by nearly 0.6m in elevation, which on flat coastal terrain translates to potentially hundreds of metres of horizontal uncertainty in where the flood edge sits. This is why engineers should always report uncertainty alongside their spatial analysis results.

• Vector contours make it easier to identify a specific elevation — you can trace the 200 m line and know exactly where that elevation occurs. The contour is a crisp, labelled feature.
• Raster DEM better shows the continuous terrain surface — the colour gradient reveals the smooth transition from valleys to ridges, showing the full 3D shape that contour lines only sample at 50 m intervals.

This illustrates the fundamental trade-off: vectors are precise for specific features, rasters capture continuous spatial variation.

What the raster preserves that vectors lose:

• Elevation values between contour lines — the raster has a value at every 8 m cell, while contours only mark specific elevation thresholds (every 50 m here)
• Local terrain detail — small ridges, depressions, or saddle points that fall between contour intervals are invisible in the vector but captured in the raster

What the vector preserves that the raster loses:

• The exact position of a known elevation — contour lines are mathematically precise boundaries, while the raster quantises position to 8 m cells
• Topological relationships — a contour line explicitly connects all points at the same elevation; in the raster, you'd need to compute this from the grid

Raster: each cell stores one elevation value. For this area: 500 x 500 = 250,000 values, one per 8 m cell. Storage is uniform — every cell costs the same regardless of terrain complexity.

Vector: each contour line stores a variable number of [x, y] vertex pairs. A complex contour in steep terrain might have hundreds of vertices; a simple one might have only a few dozen. Total vertices across all contours is much less than 250,000.

Key insight: Vectors are adaptive — they use more storage where boundaries are complex and less where they are simple. Rasters use uniform storage everywhere, even in flat areas where elevation barely changes. This means vectors are more storage-efficient for sparse features (like contours), while rasters are more efficient for truly continuous fields.

Raster DEM is the clear choice for slope calculation:

• Slope = rate of elevation change, calculated from the elevation difference between adjacent cells using a moving window (e.g., 3x3 kernel)
• The regular grid structure makes this a straightforward neighbourhood operation — compare each cell to its 8 neighbours
• Contour vectors don't have elevation values between lines, so you'd first need to interpolate them back into a grid (essentially creating a raster) before calculating slope

This is a key advantage of rasters for spatial analysis: operations like slope, aspect, flow direction, and viewshed all rely on the regular grid structure for efficient neighbourhood calculations.

Both is the ideal answer:

• DEM raster for slope/gradient analysis — compute slope everywhere to identify feasible corridors (areas below 10% gradient)
• Contour lines for route planning — engineers traditionally use contours to design switchbacks and maintain grade, as contour spacing directly shows steepness
• The DEM enables automated least-cost path analysis; contours enable manual design refinement

Resolution: For road design, 1-5 m resolution would be appropriate. The 8 m DEM here is adequate for corridor selection but too coarse for detailed design (an 8 m cell could miss a 5 m gully). Final design typically requires sub-metre LiDAR survey.

Accepting just "raster" or "both" with good justification is fine. The key insight is matching resolution to the engineering decision being made.

On the map: Greenland appears to be roughly the same size as Africa — perhaps slightly smaller.

Actual areas (students should find these online):

• Greenland: ~2,166,000 km²
• Africa: ~30,370,000 km²
• Ratio: Africa is approximately 14× larger than Greenland!

Why? Web maps use the Web Mercator projection, which preserves angles (shape) but severely distorts area, especially at high latitudes. Greenland is at ~60-84°N where the Mercator stretching factor is extreme — features are stretched by 1/cos(latitude) in each direction. At 70°N, this means features appear ~8.5× their true relative area. Africa straddles the equator where there is minimal distortion.

This is the classic "Mercator distortion" — all web maps (Google Maps, OpenStreetMap, Leaflet) inherit this problem because they use EPSG:3857 (Web Mercator).

• WGS84 (EPSG:4326) uses degrees (Longitude/Latitude). Coordinates look like (172.73°, -43.50°).
• NZTM (EPSG:2193) uses metres (Easting/Northing). Coordinates look like (1,580,000 E, 5,180,000 N).

For measuring a 100m pipe: NZTM (EPSG:2193) because it's a projected CRS with metre units. Distances can be calculated directly. In WGS84, a buffer of "100" would mean 100 degrees (~11,000 km), not 100 metres.

• set_crs() — Labels/declares the CRS of existing coordinates WITHOUT changing them. Use when data has no CRS metadata but you know what CRS the coordinates are in. "These numbers ARE in NZTM."
• to_crs() — Mathematically transforms coordinates FROM the current CRS TO a new one. Use when you need to reproject data. "Convert these NZTM coordinates TO WGS84."

Wrong use scenarios:

• Using set_crs(4326) on NZTM data — the coordinates (1575000, 5180000) get labelled as degrees. The data appears to be at longitude 1,575,000° which doesn't exist. The data is "lost" on the globe.
• Using to_crs(2193) on data that has no CRS set — GeoPandas raises an error because it doesn't know what coordinate system to transform FROM.

Students should measure and find that the same distance is expressed very differently:

• In metres: a direct, intuitive measurement (e.g., 500 m)
• In degrees: a small decimal (e.g., ~0.006°) that varies with latitude

At 43.5°S: 1° of longitude = 111,320 × cos(43.5°) ≈ 111,320 × 0.7254 ≈ 80,700 m.
So 0.001° ≈ 80.7 m.

At the equator, 0.001° ≈ 111 m. This demonstrates why degrees are unreliable for engineering measurements — the same angular distance represents different physical distances depending on latitude.

Pre-2010 datasets are now positionally inaccurate because the ground itself moved:

• Horizontal shift — features that were accurately surveyed before the earthquakes are now in the wrong location by up to several metres. Pipe networks, property boundaries, and road alignments no longer match reality.
• Vertical shift — DEMs and contour data are wrong by up to 0.5m. For flood modelling, this is critical — a 0.5m error in a flat coastal area like New Brighton means the flood boundary could be hundreds of metres off.
• Datum implications — NZGD2000 coordinates assigned before the earthquake may need transformation to match post-earthquake ground positions. LINZ released localised deformation models to handle this.

This is why temporal currency and lineage are critical data quality dimensions — the "fitness for purpose" of any pre-2010 Canterbury dataset must be reassessed.

Building counts will vary with the demo data, but typically:

• 200m buffer: a small number of buildings immediately adjacent to each marae
• 1km buffer: significantly more buildings as the circular zones overlap and cover larger areas

The buffer operation transforms a point (0-dimensional, no area) into a polygon (2-dimensional, with area). This is essential because you cannot perform area-based spatial queries (intersects, within, contains) between a point and another geometry in a meaningful way — the point has no extent.

In GeoPandas: marae.buffer(200) creates a circular polygon of radius 200m around each marae point. You can then use gpd.sjoin(buildings, marae_buffered, predicate='within') to find buildings inside those buffers.

The query geometry is a 500 m circle around the marae. The answer is about the size of each feature relative to that circle:

• Roads are long lines. A road only has to touch the circle anywhere along its length to satisfy intersects, so many roads that run past the marae match. To satisfy within the whole road must lie inside the circle, and almost no road is that short — most cross the circle's edge and continue outside. Hence a large intersects count and a small within count.
• Buildings are small polygons (tens of metres). If a building touches the circle it is almost always entirely inside it, so intersects and within agree. The only differences are the handful of buildings that straddle the circle's edge.

The lesson: the gap between intersects and within grows with feature size relative to the query geometry. For features much smaller than the buffer the two predicates are interchangeable; for features larger than, or comparable to, the buffer the choice of predicate changes the answer, so you must state which one you used.

• Euclidean distance = straight-line "as the crow flies" distance between two points
• Network distance = distance along the actual road/path network, following the routes people must travel

Example: In New Brighton, the Avon-Heathcote Estuary creates a large water body between communities. Two points on opposite sides of the estuary might be only 500m apart in Euclidean distance but require 5+ km of road travel around the estuary to actually reach each other. During a flood evacuation, using Euclidean distance would drastically underestimate response times and travel distances.

Another example: one-way streets, bridges, and road closures due to flooding mean the road network distance can be much longer than straight-line distance.

In this demo the test is feature.contains(buffer): a feature matches only if the whole marae buffer lies inside it. No building and no road is large enough to swallow a circle of 200 m to 2 km radius, so the count is always 0.

Predicate direction: contains is within with the arguments swapped. buffer.contains(building) is exactly building.within(buffer), and in the demo that returns hundreds of matches. So the predicate is not "wrong", it is being asked in the direction that can never be true.

An example from this map where contains would return matches: swap the roles so the buffer is the container: "which buildings does the 500 m buffer contain?" In GeoPandas: gpd.sjoin(marae_buffer, buildings, predicate='contains'). In general, contains is the right choice when the large polygon is on the left: suburbs containing buildings, a flood extent containing marae, a catchment containing stream reaches.

A line has length but no interior, so "within the coastline" asks whether a building sits inside something that has no inside. The question is ill-posed, not merely hard. To make it valid you first have to build a polygon from the line, for example:

• buffer the coastline by a distance ("buildings within 200 m of the coast"), or
• close the line into an area ("land" vs "sea", or a coastal hazard zone polygon), then ask which buildings are within that polygon.

Which predicates still mean something against a line:

• intersects — meaningful: does the building footprint touch or cross the line? (A road crossing the coastline would match, e.g. a jetty or causeway.)
• disjoint — meaningful: the exact opposite, no shared point.
• within — meaningless for a polygon or point feature against a line: only a point lying exactly on the line, or a piece of the line itself, can be "within" it.
• contains — meaningless here: a building would have to contain a stretch of the coastline.

The general rule: within and contains need the container to be a polygon; intersects and disjoint work between any geometry types.

Accept either "Quantile" or "Manual (flood-focused)" with sound reasoning:

• Quantile puts an equal number of features in each class. Because the elevation data is heavily skewed toward low values, the computed quantile breaks split the low coastal elevations across several classes — giving fine colour detail (and an alarming look) exactly where the land is lowest.
• Manual (flood-focused) places breaks deliberately at decision-relevant thresholds (e.g. 1m, 2m, 3m), so the low coastal strip is explicitly picked out in the highest-risk colours — arguably the most dramatic portrayal of coastal risk.

Equal Interval downplays risk because it spreads classes evenly across the full elevation range (e.g., 0-5, 5-10, 10-15, ...), lumping most of the low coastal area into a single class.

Natural Breaks (Jenks) usually gives the most balanced view, finding natural clusters in the data.

This demonstrates that classification choice is a cartographic decision that affects perception — the same data can look alarming or reassuring depending on how you classify it.

The 7 required elements:

1. Title — what the map shows
2. Legend — what colours/symbols mean
3. Scale bar — distance reference
4. North arrow — orientation
5. Data source/attribution — where the data came from
6. Projection/CRS — which coordinate system
7. Date — when the map (and its data) were produced

The demo map will typically have a title and legend but may be missing scale bar, north arrow, CRS label, and/or date. Students should identify exactly which are present and which are missing.

Rainbow (jet) colourmap is problematic because:

• Not perceptually uniform — the human eye perceives brightness changes in yellow/cyan much more than in blue/red, creating false boundaries in the data
• No intuitive order — is red "more" or "less" than green? There's no natural ordering, so viewers can't intuitively rank values
• Colour vision deficiency — approximately 8% of males have red-green colour blindness, making portions of the rainbow indistinguishable

Better choice: A sequential single-hue scheme (e.g., light blue → dark blue for flood depth) because:

• Blue intuitively represents water
• Darker = deeper is natural ordering
• Works for colour-blind viewers (varies only in lightness)

ColorBrewer "Blues" or "YlGnBu" are good options.

Raw counts are misleading because larger suburbs will always have more flooded buildings simply because they have more buildings total. A suburb with 1,000 buildings and 100 flooded (10%) looks worse than a suburb with 50 buildings and 40 flooded (80%) when mapped by count, even though the smaller suburb is proportionally much more affected.

You should map a normalised/derived variable such as:

• Percentage of buildings flooded (flooded / total × 100) — most common
• Flooded buildings per km² — accounts for different suburb areas
• Flooded building area as % of total building area — accounts for building sizes

This is the "choropleth trap" — mapping raw counts on a choropleth always favours larger areas. Always normalise.

Specific feedback:

• Colour vision issue — red-green is the most common colour blindness (deuteranopia affects ~8% of males). Roughly 1-2 councillors in the room may not be able to distinguish high risk from safe areas.
• Cultural associations — while red=danger, green=safe seems intuitive, this is a Western cultural convention. In some contexts, green can mean "action needed" and red can mean "stop/complete."
• Diverging vs sequential — red-to-green is a diverging scheme (implies a meaningful midpoint). If flood risk is a continuous scale from low to high, a sequential scheme (light to dark) is more appropriate.

Recommendation: Replace with a sequential scheme like white → orange → dark red, or use ColorBrewer "OrRd" which is colour-blind safe and perceptually ordered.

The demo estimates buffer area as total stream length × 2 × buffer distance. With ~3,241 stream reaches totalling roughly 10,000-15,000 km of channel:

• At 100m: area ≈ 12,000 km × 2 × 0.1 km ≈ ~2,400 km²
• At 500m: area ≈ 12,000 km × 2 × 0.5 km ≈ ~12,000 km²

The simple formula gives a linear relationship (5× the distance = 5× the area). However, in reality the relationship is sub-linear because:

• Buffers overlap where streams are close together — the actual buffered area is less than the sum of individual buffers (you can't count the same land twice)
• At 500m, many neighbouring stream buffers merge into a single continuous zone, so the effective new area gained per additional metre of buffer decreases
• At very large buffer distances, the riparian zone covers nearly the entire catchment and additional distance adds minimal new area

This is why the union (dissolve) of individual buffers is needed for an accurate area calculation.

• Clip — physically cuts stream geometry at the catchment boundary. Stream segments inside are kept; parts outside are removed. A stream crossing the boundary is split, and only the interior portion is retained. Use when you need the exact geometry inside the boundary (e.g., calculating total stream length within a catchment).
• Spatial join (sjoin) — tags each stream with the catchment ID it falls within, but does NOT cut the geometry. The full stream feature is kept intact. Use when you need to identify which streams belong to which catchment without modifying their geometry (e.g., listing stream IDs per catchment).

Buller example:

• "What is the total stream length inside catchment 42?" → Clip (need accurate length, streams at boundary must be cut)
• "How many stream reaches does each catchment contain?" → Spatial join (need count, not modified geometry)

When dissolving from Order 2 (3,241 polygons) to Order 3 (1,458 polygons):

• Internal boundaries disappear — where two order-2 sub-catchments shared a boundary within the same order-3 parent, that boundary is removed
• External boundaries are preserved — the outer edges of the merged catchments remain
• The result is fewer, larger polygons that tile the same landscape

Useful aggregate statistics (computed during dissolve via aggfunc):

• count — number of order-2 sub-catchments per order-3 catchment (tributary complexity)
• sum of stream length — total channel length per merged catchment
• sum of area — total catchment area (should match the merged polygon area)
• mean slope or mean elevation — average terrain characteristics per catchment

The dissolve + aggregation pattern converts detailed data into summary statistics at a coarser spatial scale.

Variation in stream count per catchment depends on:

• Catchment area — larger catchments physically contain more tributaries
• Drainage density — steep, wet areas (western ranges) have denser stream networks than flat, drier areas
• Geology — impermeable rock generates more surface drainage; permeable rock (limestone/karst) has fewer surface streams

One-to-many problem: When a stream reach intersects multiple catchment polygons (e.g., it lies exactly on a boundary), the spatial join creates duplicate rows — the same stream appears once for each catchment it touches. This inflates the stream count.

Solutions:

• Use sjoin with predicate='within' (stream must be entirely inside) rather than 'intersects'
• Use the stream's midpoint (centroid) for the join instead of the full line geometry — a point can only be within one polygon
• drop_duplicates() after the join if only one assignment is needed

1. Calculate area of each sub-catchment polygon (geometry.area in projected CRS) and add as an attribute
2. Sort catchments by area descending and select the top 10 → creates "target catchments" layer
3. Clip the stream network to the target catchments → keeps only streams inside the 10 largest catchments
4. Buffer the clipped streams by the desired riparian distance (e.g., 100m) → creates buffer polygons around each stream segment
5. Dissolve (unary_union) all buffer polygons to remove overlapping areas → single multipolygon of riparian zone
6. Calculate area of the dissolved riparian zone → total riparian area in km²

Key insight: The dissolve in step 5 is essential — without it, overlapping buffers from adjacent streams would be double-counted. The workflow uses four different vector operations (select, clip, buffer, dissolve) in sequence.

• Ridge crests: 0 (or single digits just below the crest). Nothing drains into a ridge cell. About 5% of land cells in this DEM have accumulation 0.
• Valley floors: hundreds to tens of thousands. The main stream reaching Charteris Bay collects 82,469 cells ≈ 5.3 km², the highest value in the DEM.

D8 rule: the inspector highlights the neighbour with the steepest downhill slope, i.e. the largest drop divided by distance. Diagonal neighbours are 8√2 ≈ 11.3 m away, not 8 m, so a diagonal cell needs a bigger drop to win. All the water goes to that one neighbour.

Note: zoomed in, students see long parallel arrows on smooth slopes. That is a known D8 artefact: with only 8 directions, flow on planar slopes runs in straight parallel lines instead of spreading. Multiple-flow-direction methods (D∞, MFD) address this.

There is no single right answer; a good match is roughly 500–2,000 cells (3–13 ha upstream). Reference values from the demo:

Threshold (cells)	Network length	Drainage density
100	95.5 km	6.6 km/km²
500	45.1 km	3.1 km/km²
840 (default)	35.1 km	2.4 km/km²
2,000	24.0 km	1.7 km/km²
10,000	7.6 km	0.5 km/km²

Below ~200 cells the DEM draws a "stream" in almost every gully; above ~5,000 only the main stems survive and mapped tributaries disappear. The threshold is an analyst choice, which is why it must be reported and validated.

Who is wrong? Usually neither is simply right. Many extra DEM streams are ephemeral gullies that run only in storms and that OSM volunteers never mapped; OSM is incomplete and hand-digitised. The DEM knows nothing about infiltration, springs, culverts or drains, so it can also put a channel where there is none. Validation needs field data or aerial imagery; neither layer is ground truth.

Sinks: 5 at full (float) precision, 260 at 0.1 m, 11,647 at 1 m.

The extra sinks come purely from rounding the elevations. Neighbouring cells rounded to the same value form a flat patch with no downhill neighbour, so D8 has nowhere to send the water. Integer-metre storage (as in SRTM and many older DEM products) is a classic source of artefact sinks and flats.

Catchment at the main stream's coastal outlet (with snapping):

Precision	Filled	Not filled
Float	5.28 km²	4.33 km²
0.1 m	5.36 km²	0.09 km²
1 m	5.28 km²	≈ 0

Even the 5 full-precision sinks cut the catchment by 18%, because water reaching a sink stops there and never arrives at the outlet. Filling makes the result stable at about 5.3 km² whatever the precision. The deepest fill at full precision is only 0.18 m, so these are tiny artefacts, not real ponds. Real sinks do exist (volcanic craters, quarries, karst) and should be kept, not filled.

Example on the main stream 1 km above the coast: on the channel the catchment is 4.82 km². One to five cells (8–40 m) to the side, without snapping, it is only 0.06–0.2 ha, a patch of hillslope.

A cell beside the channel sits on the valley side: its D8 arrow points into the channel, so only the strip of slope above it drains through it. All the upstream water flows past in the neighbouring channel cells without touching the pour point.

That is the "boundary through a gully" error: an outlet a few metres off the true channel captures only part of the drainage, and the boundary cuts across valleys it should contain. The fix is pour-point snapping (move the outlet to the highest-accumulation cell within a small radius, 40 m here) and always checking the result visually.

With the defaults (C = 0.35, i = 40 mm/h) at the coast: A ≈ 5.28 km², Q = 0.35 × 40 × 5.28 / 3.6 ≈ 20.5 m³/s. Exact values depend on where students click.

500 m upstream the area barely changes (≈ 5.28 km²) because no big tributary joins on that reach; by 1 km upstream it is 4.82 km² and Q ≈ 18.7 m³/s. Q scales linearly with A, so moving a structure past a tributary junction gives a step change in design flow.

Least trustworthy: generally C and i. C is a lumped guess for land cover, soil and wetness (about 0.2 for forest vs 0.9 for pavement). i depends on the return period and on a time of concentration that the rational method needs but this demo does not compute. A is the most defensible input, provided sinks were filled and the pour point snapped (Q3, Q4). The rational method is normally limited to small catchments (a few km²), so 5 km² is at its edge.

The erosion valleys radiate outward from each volcanic peak in all directions, so no single sun angle reveals all valleys equally. However:

• East (90°/30°) — best for north-south running valleys (creates cross-shadows)
• South (180°/30°) — best for east-west running valleys
• NE (45°/60°) — high sun from NE gives subtle shading on SW-facing slopes
• NW default (315°/45°) — standard cartographic convention; good general overview

Lighting direction matters because hillshade works by simulating shadows. Features aligned parallel to the light source cast no shadow and become invisible. Features perpendicular to the light show maximum contrast. This is why cartographers often create multiple hillshades or use a multidirectional hillshade for general-purpose terrain visualisation.

Steepest areas (red) are found on:

• Flanks of the volcanic peaks — the conical slopes between summit and base
• Walls of erosion valleys — where streams have cut into the volcanic flanks
• Ridge connecting the two peaks — steep sides of the saddle between summits

Flat areas (green) occur on the coastal plains to the north and the harbour basin between peaks.

Engineering applications:

• Landslide hazard mapping — slopes above ~30° are prone to failure
• Road design — maximum road gradients (typically 8-12%)
• Construction site suitability — steep terrain requires cut/fill earthworks
• Erosion risk assessment — steeper slopes = higher erosion potential (USLE equation)

Calculation:
Number of cells = 500 × 500 = 250,000 cells

Area per cell = 8m × 8m = 64 m²

Total area = 250,000 × 64 = 16,000,000 m²

Convert: 16,000,000 / 1,000,000 = 16 km²

Cross-check: the grid covers 4 km × 4 km = 16 km² ✓

• Local (cell-by-cell) — operates on each cell independently. Example: convert elevation from metres to feet (multiply each cell by 3.281). Or: reclassify elevation into zones (0-200m = lowland, 200-500m = hill, >500m = peak).
• Focal (neighbourhood) — uses a moving window of surrounding cells. Example: compute slope or hillshade using a 3×3 window — each cell's slope depends on the elevation of its 8 neighbours. Both slope and hillshade in this demo are focal operations.
• Zonal (region-based) — computes statistics within defined zones. Example: mean elevation per catchment — use catchment boundaries as zones and calculate the average DEM value within each drainage area.

(Global operations use the entire raster — e.g., computing flow accumulation where every cell's value depends on the full upstream network.)

At 100m resolution, each cell covers 10,000 m². Features lost would include:

• Narrow erosion valleys (many are only 2-4 cells wide at 8m, i.e. 16-32m) — these would be smoothed out at 100m
• Sharp ridge lines — blended with surrounding slopes
• Small-scale terrain roughness — volcanic rock outcrops, minor gullies
• Peak summits — the highest point would be averaged with surrounding lower cells, reducing the apparent peak height

Resampling methods:

• DEM (continuous data) → Bilinear or bicubic interpolation — produces smooth interpolated values. Appropriate because elevation varies continuously and intermediate values are physically meaningful (a value of 450m between a 400m and 500m cell makes sense).
• Land use (categorical data) → Nearest neighbour — assigns each new cell the value of the closest original cell. Appropriate because land-use categories are discrete — there is no meaningful "average" of "forest" and "urban." Bilinear would create nonsensical intermediate values.

At 50m threshold: the demo shows the count of cells with elevation > 0 and < 50m. The exact number depends on the DEM data, but typically ~20,000–30,000 cells.

Area calculation:
Area = cell_count × 8² = cell_count × 64 m²
e.g. if 25,000 cells: Area = 25,000 × 64 = 1,600,000 m² = 1.60 km²

At 100m threshold: the flooded area does NOT simply double. This is because:

• Elevation distribution is not uniform — the hypsometric curve of Banks Peninsula is nonlinear
• The terrain has steep volcanic flanks between ~50m and ~200m, meaning relatively fewer cells in that elevation band
• Most low-elevation cells are on the coastal plains, which are already captured by the 50m threshold
• This is a local raster operation — each cell is independently classified as flooded or not based solely on its own elevation value

The relationship between threshold and flooded area is nonlinear and depends on the terrain's hypsometric distribution (the cumulative distribution of area with elevation).

In this demo the three components are formulas computed per elevation-contour line and combined into one colour-coded score — they are not separate toggleable layers:

• Hazard = modelled flood depth: contour elevation vs the selected sea level. Visible as the dashed flood boundary — higher sea level = deeper flooding = higher hazard.
• Exposure = an elevation-based proxy for the built environment (low-lying flat land carries most buildings and roads). The roads and coastline layers on the map are the visible context for this — note the densest road network sits on the lowest ground.
• Vulnerability = proximity to the coast combined with low elevation — assets nearer the shoreline score as more susceptible.

Beyond the demo: a full MCE would use separate data layers for each component — flood depth rasters (hazard); building footprints, road networks, marae and infrastructure (exposure); building age/construction type, fragility curves, or NZDep social deprivation (vulnerability). Note that buildings and marae layers are NOT in this demo — students should recognise the demo substitutes elevation proxies for real asset data.

Risk is spatial — all three components vary across the map, so the risk score varies by location. GIS is the platform that combines them.

With Hazard=60%: the risk map is dominated by flood depth — areas with deeper flooding show high risk regardless of how many buildings are there. The map highlights WHERE the water goes.

With equal weights: risk is more balanced — an area with moderate flooding but high building density and vulnerable populations may score higher than a deeply flooded but empty area. The map highlights WHERE the consequences are.

For coastal flood planning: Hazard-weighted (60/30/10) is more appropriate for emergency response planning (need to know where water goes first). Equal weights are more appropriate for long-term land use planning (need to balance all factors).

There is no single "correct" weighting — it depends on the decision being made. This is why MCE always requires sensitivity analysis.

At 2m sea level rise (read directly from the demo's impact panel):

• Flooded area: 2.39 km²
• Flooded houses: 2.39 × 150 = 359 houses (panel shows 359)
• Flooded roads: 2.39 × 10 = 23.9 km
• House loss: 359 × $600,000 = ~$215.4 million NZD
• Road loss: 23.9 × $1,500,000 = ~$35.9 million NZD
• Total: ~$251.3 million NZD

Are assumptions reasonable?

• $600k/house is low for 2025 Christchurch — median house price is closer to $650-750k. However, flood damage rarely destroys 100% of value, so a lower figure may represent partial damage.
• $1.5M/km of road is reasonable for urban road replacement including services underneath.
• Missing: contents insurance, business disruption, infrastructure (pipes, power), health costs, displacement costs. Actual losses could be 2-3× higher.

Cost comparison (demo mitigation cards — relocation covers ALL flooded houses):

• Seawall: 12 km × $50M/km = $600 million NZD at 1-2m (the +$10M/km per metre surcharge kicks in above 2m: $720M at 3m, $840M at 4m)
• Managed retreat at 2m: 359 flooded houses × 4 people × $264,450/person ≈ $380 million NZD
• At other levels: ~$9.5M at 1m, ~$529M at 3m — in this demo retreat is cheaper than the wall at every scenario

Relocation looks cheaper, but the comparison isn't straightforward:

• Seawall protects ALL assets behind it (thousands of buildings, roads, infrastructure) — cost per protected asset is actually lower
• Seawall has ongoing maintenance costs and may fail catastrophically under extreme events
• Relocation is permanent — removes exposure entirely, with no residual risk of defence failure
• Relocation has social/cultural costs that are hard to quantify (community disruption, tangata whenua connections to land)

Extension discussion (not shown in the demo): a targeted policy might relocate only the ~100 most vulnerable houses: 100 × 4 × $264,450 ≈ $105.8M — far cheaper than either option above, but it leaves the remaining exposure unaddressed.

Recommendation: A hybrid approach — relocate the most at-risk houses in the lowest-lying areas AND invest in targeted coastal protection for the remaining community. Pure cost-benefit analysis is insufficient; social and cultural factors must be weighed.

Sensitivity analysis tests how much the output changes when inputs vary. It's essential because:

• MCE weights are subjective — different stakeholders would assign different weights
• If the result changes dramatically with a small weight adjustment, the conclusion is not robust and shouldn't be used for major decisions without further analysis
• If the result is stable across a range of weights, the conclusion is robust and can be used with more confidence

In the demo: If changing Hazard weight from 50% to 30% causes the #1 risk area to become #3, the ranking is sensitive to hazard weighting. This means the team needs to invest more effort in agreeing on hazard weight, or acknowledge the uncertainty in their recommendation.

Best practice: Report results as ranges ("Area X ranks #1-3 depending on weight scenario") rather than single values, and use equal-weight baseline as a reference.

Note: in this demo all three criteria derive from elevation and distance-to-coast, so they are strongly correlated — students should EXPECT the highest-risk ranking to be robust (barely move) under ±20% weight changes. That is itself the teaching point: correlated criteria mute sensitivity, whereas independent criteria (e.g. real building-stock exposure) would swing more.