Cheat Sheet – Chapter III, Exercise 1.1

This document summarises the key concepts and coding patterns from Exercise 1.1. Code snippets are shown for illustration only, so they are not executed.

1. What is a Web Coverage Service (WCS)?

A Web Coverage Service (WCS) is an open standard from the Open Geospatial Consortium (OGC) that allows you to access raw raster data (such as elevation models or paleogeographic maps) directly over the web, without downloading files manually.

Unlike a Web Map Service (WMS), which returns an image, a WCS returns the actual numeric values of the raster, so you can do real computations on the data.

Anatomy of a WCS URL

A WCS request is built from a set of key-value parameters joined into a URL. Every parameter has a specific role:

Parameter	Role	Example
`base_url`	Root address of the server	`https://geoserver.panalesis.org/geoserver/`
`workspace`	Data store on the server	`panalesis_atlas`
`service`	Which OGC service to use	`WCS`
`version`	Protocol version	`1.0.0`
`request`	Type of request	`GetCoverage`
`coverage`	Layer identifier	`workspace:layer_name`
`crs`	Coordinate Reference System	`EPSG:4326` or `EPSG:54034`
`bbox`	Spatial extent (minX, minY, maxX, maxY)	`-180,-90,180,90`
`resx` / `resy`	Pixel size in CRS units	`0.1` (degrees) or `10000` (metres)
`time`	Time step (ISO 8601)	`2100-01-01T00:00:00.000Z`
`format`	Output format	`GEOTIFF`

Constructing the URL in Python

base_url      = "https://geoserver.panalesis.org/geoserver/"
workspace     = "panalesis_atlas"
service_type  = "WCS"
service_version = "1.0.0"
request_type  = "GetCoverage"
layer_name    = "palaeogeography"
crs           = "EPSG:4326"
bbox          = "-180,-90,180,90"
resx          = "0.1"
resy          = "0.1"
time          = "2100-01-01T00:00:00.000Z"   # encodes 100 Ma
format        = "GEOTIFF"

wcs_url = (
    rf"{base_url}{workspace}/{service_type.lower()}?"
    rf"service={service_type}&"
    rf"version={service_version}&"
    rf"request={request_type}&"
    rf"coverage={workspace}:{layer_name}&"
    rf"crs={crs}&"
    rf"bbox={bbox}&"
    rf"resx={resx}&"
    rf"resy={resy}&"
    rf"time={time}&"
    rf"format={format}"
)

Time encoding trick

Geospatial data servers (like GeoServer) follows ISO 8601. Geological ages in Ma are encoded by adding the age to the year 2000. For example, 100 Ma → year 2100 → 2100-01-01T00:00:00.000Z. Present-day (0 Ma) → 2000-01-01T00:00:00.000Z.

2. Loading a WCS Layer and Its Metadata

When you load a layer from a WCS, you receive not only pixel values you also get a rich set of metadata describing the spatial properties of the raster.

import requests
import numpy as np
from rasterio.io import MemoryFile

data = requests.get(wcs_url).content  # raw bytes from the server

with MemoryFile(data) as memfile:
    with memfile.open() as dataset:

        # --- Read the pixel values ---
        raster = dataset.read(1).astype(float)

        # --- Read the metadata ---
        print(dataset.meta)       # driver, dtype, nodata, width, height, CRS, transform
        print(dataset.bounds)     # spatial extent: left, bottom, right, top
        print(dataset.crs)        # coordinate reference system
        print(dataset.transform)  # affine transformation matrix

Key metadata attributes

Attribute	What it tells you
`meta`	Full summary: format, data type, no-data value, pixel dimensions, CRS, transform
`bounds`	Spatial extent of the raster (in CRS units: degrees or metres)
`crs`	The coordinate reference system (e.g. EPSG:4326)
`transform`	Affine matrix mapping pixel indices (row, col) → real-world coordinates
`nodata`	Sentinel value for missing data — must be replaced with `np.nan` before analysis

The `transform` attribute in detail

The affine transform is a 3×3 matrix that converts pixel positions to geographic coordinates. It encodes:

The pixel size (resolution) in x and y
The origin of the raster (top-left corner coordinates)
Any rotation (usually zero for north-up rasters)

# Example transform output:
# | pixel_width   0          x_origin |
# | 0            -pixel_height y_origin |
# | 0             0           1        |

# Accessing pixel size:
pixel_width  = dataset.transform[0]   # e.g.  0.1 degrees (or 10000 m)
pixel_height = abs(dataset.transform[4])  # e.g. 0.1 degrees (note: negative in rasters)

Why does this matter?

Pixel area = pixel_width x pixel_height. You need the transform to convert a pixel count into a real-world area.

3. What Makes a Good Map?

A raw imshow() call is not a map. A proper map requires four elements:

Element	Why it matters	How to add it
Coherent symbology	Colours should be perceptually meaningful (e.g. blue = ocean, green/brown = land)	Custom `LinearSegmentedColormap` + `Normalize`
Axis labels with units	Readers must know what the axes represent	`ax.set_xlabel()` / `ax.set_ylabel()`
Legend / colour bar	Links colours to values	`plt.colorbar(img, ax=ax)` + `cbar.set_label()`
Contextual information	Title, date, data source — the audience needs to know what they are looking at	`ax.set_title()`

Minimal example of a well-formed map

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

# 1. Define a meaningful colour ramp
colormap_entries = [
    (-9000, "#001e64"),
    (-4000, "#136ca0"),
    (    0, "#006147"),
    (  200, "#107b30"),
    ( 1000, "#e8d67d"),
    ( 3000, "#821e1e"),
    ( 9000, "#ffffff"),
]
elevations = [e[0] for e in colormap_entries]
colors     = [e[1] for e in colormap_entries]
vmin, vmax = elevations[0], elevations[-1]
norm_vals  = [(e - vmin) / (vmax - vmin) for e in elevations]
cmap = mcolors.LinearSegmentedColormap.from_list(
    "palaeogeography", list(zip(norm_vals, colors))
)
norm = mcolors.Normalize(vmin=vmin, vmax=vmax)

# 2. Plot with extent so axes show real coordinates
fig, ax = plt.subplots()
img = ax.imshow(raster, cmap=cmap, norm=norm, extent=extent, origin="upper")

# 3. Add all required map elements
cbar = plt.colorbar(img, ax=ax)
cbar.set_label("Elevation (m)")          # legend label
ax.set_xlabel("Longitude (°)")           # x-axis label
ax.set_ylabel("Latitude (°)")            # y-axis label
ax.set_title("Palaeogeographic map – 100 Ma (PANALESIS)")  # context

plt.tight_layout()
plt.show()

4. Coordinate Reference Systems and Area Calculations

Never compute areas in EPSG:4326 (lat/lon)

When a map uses geographic coordinates (degrees), the pixel size is not constant in meters: a 0.1x0.1° pixel near the equator is much larger than the same pixel near the poles. Area calculations based on pixel counts will be wrong.

Solution: use an equal-area projection

Switch to a CRS that preserves area (e.g.ESRI:54034: World Cylindrical Equal Area). In this CRS, all pixel sizes are equal in square metres, so multiplying pixel count by pixel area gives a correct result.

# Equal-area WCS request
crs  = "EPSG:54034"
bbox = "-20037508.34,-6363885.33,20037508.34,6363885.33"  # in metres
resx = "10000"  # 10 km pixels
resy = "10000"

CRS	Units	Use for
EPSG:4326	Decimal degrees	Displaying maps, basic visualisation
EPSG:54034	Metres (equal area)	Computing areas, statistics

5. Writing a Function in Python

A function bundles a reusable piece of logic. Every function has three parts:

Arguments: the inputs it needs
Body: what it does
Return value: what it gives back when called

def get_land_ocean_area(data, transform, extent, plot=False):
    """
    Compute land and ocean areas from an elevation raster.

    Arguments
    ---------
    data: 2-D numpy array of elevation values (in metres)
    transform: rasterio affine transform of the dataset
    extent: [xmin, xmax, ymin, ymax] used only when plot=True
    plot: bool, whether to display the land/ocean mask

    Returns
    -------
    land_area: float, total land area in m²
    ocean_area: float, total ocean area in m²
    """
    # Pixel area from the transform (units², here m²)
    pixel_area = abs(transform[0] * transform[4])

    # Boolean masks
    land_mask  = data >= 0    # True where land
    ocean_mask = data <  0    # True where ocean

    if plot:
        import matplotlib.pyplot as plt
        fig, ax = plt.subplots()
        ax.imshow(land_mask, cmap="Greys", extent=extent, origin="upper")
        ax.set_xlabel("Easting (m)")
        ax.set_ylabel("Northing (m)")
        ax.set_title("Land / Ocean mask")
        plt.tight_layout()
        plt.show()

    # Area = number of True pixels × area of one pixel
    land_area= float(np.sum(land_mask)  * pixel_area)
    ocean_area= float(np.sum(ocean_mask) * pixel_area)

    return land_area, ocean_area   # <-- two values returned as a tuple

Calling the function

land_area, ocean_area = get_land_ocean_area(
    data= raster,
    transform= dataset.transform,
    extent= extent,
    plot= True
)
print(f"Land area= {land_area} m²")
print(f"Ocean area= {ocean_area} m²")

Good function design habits

Give descriptive argument names
Include a short text at the beginning of the funciton explaining inputs and outputs
Keep functions focused on one task

6. `GetCoverage`: How it Works

GetCoverage is the WCS request type that downloads raster data. There are two other common request types:

Request type	What it does
`GetCapabilities`	Returns a description of what the server offers (layers, formats, CRS…)
`DescribeCoverage`	Returns metadata for a specific layer (spatial/temporal extent, resolution…)
`GetCoverage`	Downloads the actual raster data for a given bbox, CRS, and time

In Python, the response is raw bytes, which rasterio.MemoryFile can open without writing to disk:

from rasterio.io import MemoryFile

response_bytes = requests.get(wcs_url).content   # download

with MemoryFile(response_bytes) as memfile:       # wrap in memory
    with memfile.open() as dataset:               # open as rasterio dataset
        raster = dataset.read(1).astype(float)    # read band 1
        transform = dataset.transform
        bounds = dataset.bounds

7. From a Single Time Step to a `for` Loop

Once your code works for a single map, it is very easy to repeat it for all available time steps using a for loop.

Step 1: Discover available time steps

from xml.etree import ElementTree as ET

def get_available_times(base_url, workspace, service_type, service_version, layer_name):
    """Query the server for all available time steps."""
    describe_url = (
        f"{base_url}{workspace}/{service_type.lower()}?"
        f"service={service_type}&version={service_version}&"
        f"request=DescribeCoverage&coverage={workspace}:{layer_name}"
    )
    response = requests.get(describe_url)
    root = ET.fromstring(response.content)   # parse the XML response

    # Extract all <timePosition> values
    times = sorted(set(
        el.text for el in root.iter()
        if el.tag.endswith('timePosition') and el.text
    ))
    return sorted(times)

times = get_available_times(base_url, workspace, service_type, service_version, layer_name)

Step 2: Loop over time steps and accumulate results

results = {}   # dictionary to store results keyed by age

for time in times:
    # Decode geological age from the ISO timestamp
    geological_age = int(time[:4]) - 2000  # e.g. "2100-..." → 100 Ma

    # Build the WCS URL for this time step
    wcs_url = (
        rf"{base_url}{workspace}/{service_type.lower()}?"
        rf"service={service_type}&version={service_version}&"
        rf"request={request_type}&coverage={workspace}:{layer_name}&"
        rf"crs={crs}&bbox={bbox}&resx={resx}&resy={resy}&"
        rf"time={time}&format={format}"
    )

    # Download and analyse
    data = requests.get(wcs_url).content
    with MemoryFile(data) as memfile:
        with memfile.open() as dataset:
            raster = dataset.read(1).astype(float)
            nodata = dataset.nodata
            if nodata is not None:
                raster[raster == nodata] = np.nan
            bounds = dataset.bounds
            extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]

    land_area, ocean_area = get_land_ocean_area(raster, dataset.transform, extent)

    # Store results in the dictionary
    results[geological_age] = {
        'land_area':  land_area,
        'ocean_area': ocean_area
    }

    del raster, data   # free memory before next iteration

8. Storing Results in a Dictionary

A dictionary (dict) maps keys to values. It is ideal for accumulating results inside a loop when each iteration produces one record.

# Empty dictionary — created before the loop
results = {}

# Inside the loop, add an entry with the key being the geological age, and the values being the land and ocean areas.
results[geological_age] = {
    'land_area':  land_area,
    'ocean_area': ocean_area
}

# After the loop, access values like this:
print(results[100])            # → {'land_area': ..., 'ocean_area': ...}
print(results[100]['land_area'])  # → numeric value

Extracting lists for plotting

ages = list(results.keys())
land_areas = [results[age]['land_area']  for age in ages]
ocean_areas = [results[age]['ocean_area'] for age in ages]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

ax1.plot(ages, land_areas,  color='darkgreen', marker='o', markersize=4)
ax1.set_xlabel('Age (Ma)')
ax1.set_ylabel('Land area (m²)')
ax1.set_title('Land area through the Phanerozoic')
ax1.invert_xaxis()   # geological convention: oldest on the left

ax2.plot(ages, ocean_areas, color='darkblue',  marker='o', markersize=4)
ax2.set_xlabel('Age (Ma)')
ax2.set_ylabel('Ocean area (m²)')
ax2.set_title('Ocean area through the Phanerozoic')
ax2.invert_xaxis()

plt.tight_layout()
plt.show()

invert_xaxis()

In Earth sciences, time axes conventionally run from oldest on the left to youngest (present-day) on the right. ax.invert_xaxis() enforces this.

Summary Cheatsheet

Concept	Key idea	Code pattern
WCS URL	Assemble parameters into a query string	`f"...?service=WCS&request=GetCoverage&..."`
Load raster	`requests` → bytes → `MemoryFile` → `dataset`	`MemoryFile(data) as mf: mf.open() as ds`
Metadata	`bounds`, `crs`, `transform`, `nodata`	`dataset.meta`, `dataset.transform`
Good map	Cmap + norm + extent + labels + title	`imshow(..., extent=...) + colorbar + labels`
Area in m²	Use equal-area CRS; pixel count × pixel area	`EPSG:54034`; `abs(t[0]t[4]) np.sum(mask)`
Function	Args → body → `return`	`def f(x, y, plot=False): ... return a, b`
Loop	Repeat for all time steps	`for time in times: ...`
Dictionary	Store keyed results	`results[age] = {'land': ..., 'ocean': ...}`