How did Earth look 100 million years ago, and how has its geography changed over hundreds of millions of years?
In this module, you will:
Start simple: Explore a single paleogeographic map and visualize Earth at 100 Ma.
Compare worlds: Compute and interpret differences in land and ocean distribution between two time slices.
Go deep time: Extend your analysis to the entire Phanerozoic (last 540 million years) and identify major trends.
Earth’s geography controls climate, ocean circulation, and biodiversity. By reconstructing these changes, you’ll learn how scientists use data and coding to answer fundamental questions about our planet’s history.
Find the data: Search on the web for maps of the Earth’s past.
Load and inspect: Read raster data in Python, check metadata, and understand its structure.
Compute basic stats:
3.1. Land vs. ocean fraction for a chosen time slice.
3.2. Differences between two slices.
3.3. A time series of land fraction across the Phanerozoic.
Visualize:
4.1 Global maps using matplotlib and cartopy.
4.2. A trend plot showing how land/ocean distribution evolved.
transform attribute. What does this attribute describe ?This tutorial requires the following libraries:
requestsmatplotlibrasterionumpyYou can check if they are installed (and check which version you have) by running:
libs = ["requests", "matplotlib", "rasterio", "numpy"]
for lib in libs:
try:
module = __import__(lib)
print(lib + " is installed")
except ImportError:
print(lib + " is NOT installed")requests is installed
matplotlib is installed
rasterio is installed
numpy is installedIf you encounter an error, come back to Chapter II: Setting up your Python environment and follow the instructions.
You want to learn what the Earth looked like 100 millions of years ago. Where do you look for such kind of data ?
Using a search engine, results will typically yield something like:
→ Have a look at the various results, how many are actual maps (e.g. with formats that are compatible with geographic information systems (GIS) software ?
Do some maps look similar ? Do they share a common course ? Or, on the contrary, do they look drastically different ?
In our exercise, we are particularly interested in having maps of the entire Earth, and we we will see further, not only one (unique) map but a series of maps showing us the evolution of the Earth for hundreds of millions of years.
Among the most well known sources of data, you have a dataset created by C. Scotese & N. Wright as part of the PALEOMAP project, accessible here https://doi.org/10.5281/zenodo.5460860
This dataset contains 88 maps showing the palaeo Digital Elevation Model (palaeo-DEM) over the Phanerozoic (from 541 Ma to present-day).
For 100Ma, the map would look like (with map resolution at 1x1°):
→ Try and display the same map in your Desktop GIS.
→ Explore the map. → Particularly, pay attention to the rendering of the topography. What do you observe on continents ? What about oceans ?
It is easier to get the maps at 1x1° resolution from the EarthByte webpage listed in the Zenodo webpage record.
We found our first example very useful to get a sense of how the Earth looked like 100 millions of years ago. However, we noticed some limitations. We also want to benefit from the experience we have at the University of Geneva in this domain.
We are going to use the PANALESIS plate tectonic model and its derived palaeogeographic maps, that were developed at UNIGE.
The data is also available on Zenodo at https://doi.org/10.5281/zenodo.15396265
→ Load the 100Ma map and compare it to the Scotese & Wright map. The PANALESIS map should look something like:
In the two examples above, we manually downloaded maps from an open data archive, and opened a single one with a dekstop GIS software. This may be useful to have a first look at some data, yet, if we want to get a deeper understanding of what our data can reveal, we need a simpler (more direct) way to access and most importantly, to analyze our data.
Fortunately for us, there are some easy ways to serve geospatial data directly to users via the web, to avoid downloading data on your local machine each time you need it. To do this, we will have a look at the Open Geospatial Consortium (OGC) web services, that include the Web Map Servives (WMS), or the Web Coverage Services (WCS). These two types of services are used on raster data. For vector data, the dedicated service is the Web Feature Services (WFS).
We will use the OGC web services to get data that is stored in a server dedicated to serve geospatial data (GeoServer). The common way to get data is through a request using an URL that contains all elements required by the web services. In our case, we will need to specify some information on the following elements:
| Parameter | Description | Example value | Required |
|---|---|---|---|
base_url |
Root URL of the GeoServer instance hosting the data | https://geoserver.panalesis.org/geoserver/ |
Yes |
workspace |
Name of the data store (workspace) inside the server | panalesis_atlas_epsg_4326 |
Yes |
service_type |
OGC service to use — WMS or WCS for raster data, WFS for vector data | WCS |
Yes |
service_version |
Version of the OGC service protocol | 1.0.0 |
Yes |
request_type |
Type of request sent to the service (e.g. GetCoverage, GetMap, GetFeature) | GetCoverage |
Yes |
layer_name |
Name of the specific layer to access within the workspace | palaeogeography_4326 |
Yes |
crs |
Coordinate Reference System — defines the spatial projection | EPSG:4326 |
Yes |
bbox |
Bounding box as minLon, minLat, maxLon, maxLat — defines the spatial extent | -180,-90,180,90 |
Yes |
resx / resy |
Spatial resolution in CRS units (degrees for EPSG:4326) | 0.1 |
Yes |
time |
ISO 8601 timestamp to select a specific time step from a temporal dataset | 2100-01-01T00:00:00.000Z |
Optional |
format |
Desired output format for the downloaded data | GEOTIFF |
Yes |
We are going to retrieve the data from the PANALESIS server, which is hosted by the University of Geneva, and available at https://geoserver.panalesis.org/geoserver/web/?2.
If you open this link in your browser, and click on “Layer Preview”, you will see all available layers. We will try and explore (almost) all of them in this course. For now, our goal is still to access the palaeogeographic map at 100Ma.
Typically, in our example, a correct URL for a web coverage service (WCS), which will return a full map will look like:
base_url = "https://geoserver.panalesis.org/geoserver/"
workspace = "panalesis_atlas_epsg_4326"
service_type = "WCS"
service_version = "1.0.0"
request_type = "GetCoverage"
layer_name = "palaeogeography_4326"
crs = "EPSG:4326"
bbox = "-180,-90,180,90"
resx = "0.1"
resy = "0.1"
time = "2100-01-01T00:00:00.000Z"
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}"
)
print(wcs_url)https://geoserver.panalesis.org/geoserver/panalesis_atlas_epsg_4326/wcs?service=WCS&version=1.0.0&request=GetCoverage&coverage=panalesis_atlas_epsg_4326:palaeogeography_4326&crs=EPSG:4326&bbox=-180,-90,180,90&resx=0.1&resy=0.1&time=2100-01-01T00:00:00.000Z&format=GEOTIFF→ Test what happens when you copy and paste this link in your browser.
Now we would like to directly retrieve the data, we can simply associate it to a temporary layer and visualize it.
import requests
import matplotlib.pyplot as plt
from rasterio.io import MemoryFile
data = requests.get(wcs_url).content
with MemoryFile(data) as memfile:
with memfile.open() as dataset:
plt.imshow(
dataset.read(1),
cmap='terrain'
)
plt.colorbar()
plt.show()
You can now recognize the same kind of map as shown above, with two major differences. What are they ?
We are going to resolve this by adding a correct color ramp and setting the axes to reflect longitude and latitude instead of pixel index.
import matplotlib.colors as mcolors
import numpy as np
colormap_entries = [
(9000, "#ffffff", "9000"),
(7000, "#cecece", "7000"),
(5000, "#a1a1a1", "5000"),
(3000, "#821e1e", "3000"),
(2000, "#a34400", "2000"),
(1000, "#e8d67d", "1000"),
(200, "#107b30", "200"),
(0, "#006147", "0"),
(-100, "#b0e2ff", "-100"),
(-500, "#87cefa", "-500"),
(-2000, "#188ccd", "-2000"),
(-4000, "#136ca0", "-4000"),
(-7000, "#003266", "-7000"),
(-9000, "#001e64", "-9000"),
(-11000, "#000050", "-11000"),
]
colormap_entries_sorted = sorted(
colormap_entries,
key=lambda x: x[0]
)
elevations = [e[0] for e in colormap_entries_sorted]
colors = [e[1] for e in colormap_entries_sorted]
vmin, vmax = elevations[0], elevations[-1]
norm_values = [(e - vmin) / (vmax - vmin) for e in elevations]
cmap = mcolors.LinearSegmentedColormap.from_list(
"palaeogeography",
list(zip(norm_values, colors))
)
norm = mcolors.Normalize(vmin=vmin, vmax=vmax)
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]
fig, ax = plt.subplots()
img = ax.imshow(
raster,
cmap=cmap,
norm=norm,
extent=extent,
origin="upper"
)
cbar = plt.colorbar(img, ax=ax)
cbar.set_label("Elevation (m)")
ax.set_xlabel("Longitude (°)")
ax.set_ylabel("Latitude (°)")
plt.show()
We now have the same result as the one above, and all elements making it a comprehensive map are present (realistic colors, legend, axes labels). However, we can go further and read some basic metadata from the WCS request itself.
data = requests.get(wcs_url).content
with MemoryFile(data) as memfile:
with memfile.open() as dataset:
print(dataset.meta)
print(dataset.bounds)
print(dataset.crs)
print(dataset.transform){'driver': 'GTiff', 'dtype': 'float32', 'nodata': -9999.0, 'width': 3600, 'height': 1800, 'count': 1, 'crs': CRS.from_wkt('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'), 'transform': Affine(0.1, 0.0, -180.0,
0.0, -0.1, 90.0)}
BoundingBox(left=-180.0, bottom=-90.0, right=180.0, top=90.0)
EPSG:4326
| 0.10, 0.00,-180.00|
| 0.00,-0.10, 90.00|
| 0.00, 0.00, 1.00|These few lines will give you useful information to understand the nature of the data you are dealing with. the result will inform you about the full list of attributes, including the data format, type, how no-data value pixels are represented, the raster dimensions, the coordinate reference system, and the “transform” which is the affine transformation matrix to match the pixel position (x and y) in the original file with their corresponding longitude and latitude.
→ Before going to the next part, write a quick summary about the data we just discuss:
As you might have observed, we can immediately spot differences between the world at 100Ma compared to the world we know today. We however need a way to quantify these differences objectively. Fortunately for us, this ca be performed quite easily with a few functions that we will explore now.
Before diving into the analysis of our map, we however need to recognize a limit with its dimensions. As noted before, the current CRS are the “traditional” EPSG:4326 (WGS 1984), with units expressed in decimal degrees. This is not useful if we want to compute areas. For this reason, we will switch to the same map that are represented using an equal area projection. In our case, we will use an equal area projection.
The new WCS url is defined as follows:
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:54034"
bbox = "-20037508.34,-6363885.33,20037508.34,6363885.33"
resx = "10000"
resy = "10000"
time = "2100-01-01T00:00:00.000Z"
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}"
)
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]
fig, ax = plt.subplots()
img = ax.imshow(
raster,
cmap=cmap,
norm=norm,
extent=extent,
origin="upper"
)
cbar = plt.colorbar(img, ax=ax)
cbar.set_label("Elevation (m)")
ax.set_xlabel("Easting (m)")
ax.set_ylabel("Northing (m)")
plt.show()
Notice how the bounds value changed: everything is not expressed in meters, which are much better units than decimal degrees when calculating areas.In our example, we will have a quick look at the land and oceanic area, which are simply defined as the count of pixels respectively above, and below or equal to 0m of elevation, multiplied by the area of pixels.
This can be implemented using a very simple function:
def get_land_ocean_area(data, transform, extent, plot=False):
pixel_area = abs(transform[0] * transform[4])
# Masks
land_mask = data >= 0
ocean_mask = data < 0
if plot:
fig, ax = plt.subplots()
img = ax.imshow(
land_mask,
cmap="Greys",
interpolation="none",
extent=extent,
origin="upper"
)
# no colorbar for a mask
ax.set_xlabel("Longitude (m)")
ax.set_ylabel("Latitude (m)")
plt.tight_layout()
plt.show()
land_area = np.sum(land_mask) * pixel_area
ocean_area = np.sum(ocean_mask) * pixel_area
return land_area, ocean_areaThis function will take the input raster data and simply look for pixels above sea-level to define them as land, and do the same with pixels at or below sea-level and define in the ocean mask. Then, it will calculate the respective areas by multiplying them by the pixel area.
You can call this function on the raster we previously loaded:
land_area, ocean_area = get_land_ocean_area(
raster,
dataset.transform,
extent,
plot=True
)
print("Land area =", land_area, "m²")
print("Ocean area =", ocean_area, "m²")
Land area = 135611979059946.42 m²
Ocean area = 374453636909075.7 m²→ What do these number mean ? Are they realistic ?
We want to compare this value to the present-day world. Now, we know how to do it, we simply need to change the data source and set it to the present-day reconstruction. You might have observed how the time dimension was defined in the previous requests. For instance, the 100Ma map was obtained using time = "2100-01-01T00:00:00.000Z". This is a trick to follow the ISO 8601 standard. We simply add the geological age (in millions of years) to the year 2000. For the 100Ma map, the code is therefore 2100-01-01T00:00:00.000Z.
To obtain the map at present-day, we can create a new WCS request with an updated time value and calculate the same statistics as we did for 100Ma:
time = "2000-01-01T00:00:00.000Z"
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}"
)
data = requests.get(wcs_url).content
with MemoryFile(data) as memfile:
with memfile.open() as dataset:
raster = dataset.read(1).astype(float)
land_area, ocean_area = get_land_ocean_area(
raster,
dataset.transform,
extent,
plot=True
)
print("Land area =", land_area, "m²")
print("Ocean area =", ocean_area, "m²")
Land area = 158193714972089.66 m²
Ocean area = 351871900996932.5 m²→ Are these number realistic ? Search for the actual land and ocean areas and compare them with your results.
Finally, we are going to extend this analysis to all our available maps. For this, we will simply create a loop that iterates through all time steps and extract the values for land and ocean areas. First, we can define a function that lists all available reconstructions. We can do this using the The xml.etree.ElementTree module, a Python standard library tool for parsing and navigating XML documents as a tree structure, where each element (tag, attributes, text) is accessible as a node.
from xml.etree import ElementTree as ET
def get_available_times(base_url, workspace,
service_type, service_version,
layer_name):
describe_url = (
f"{base_url}{workspace}/{service_type.lower()}?"
f"service={service_type}&"
f"version={service_version}&"
f"request=DescribeCoverage&"
f"coverage={workspace}:{layer_name}"
)
response = requests.get(describe_url)
root = ET.fromstring(response.content)
times = sorted(set(
el.text for el in root.iter()
if el.tag.endswith('timePosition') and el.text
))
return sorted(times)The response from the request we send to the server comes back as XML, a nested syntax with tags that is very useful for machines to read text, but not so much for us. In our case, we sent a DescribeCoverage for a given layer, which yields back information about spatial and temporal coverage of the layer. We are only interested in the time property so this is why we filter out only values inserted inside the timePosition tag.
If we call the function:
times = get_available_times(
base_url,
workspace,
service_type,
service_version,
layer_name
)
print(times)['2000-01-01T00:00:00.000Z', '2006-01-01T00:00:00.000Z', '2011-01-01T00:00:00.000Z', '2015-01-01T00:00:00.000Z', '2020-01-01T00:00:00.000Z', '2033-01-01T00:00:00.000Z', '2040-01-01T00:00:00.000Z', '2048-01-01T00:00:00.000Z', '2056-01-01T00:00:00.000Z', '2068-01-01T00:00:00.000Z', '2084-01-01T00:00:00.000Z', '2094-01-01T00:00:00.000Z', '2100-01-01T00:00:00.000Z', '2113-01-01T00:00:00.000Z', '2120-01-01T00:00:00.000Z', '2133-01-01T00:00:00.000Z', '2140-01-01T00:00:00.000Z', '2154-01-01T00:00:00.000Z', '2165-01-01T00:00:00.000Z', '2180-01-01T00:00:00.000Z', '2200-01-01T00:00:00.000Z', '2210-01-01T00:00:00.000Z', '2220-01-01T00:00:00.000Z', '2230-01-01T00:00:00.000Z', '2240-01-01T00:00:00.000Z', '2250-01-01T00:00:00.000Z', '2270-01-01T00:00:00.000Z', '2290-01-01T00:00:00.000Z', '2300-01-01T00:00:00.000Z', '2315-01-01T00:00:00.000Z', '2331-01-01T00:00:00.000Z', '2350-01-01T00:00:00.000Z', '2370-01-01T00:00:00.000Z', '2383-01-01T00:00:00.000Z', '2393-01-01T00:00:00.000Z', '2408-01-01T00:00:00.000Z', '2420-01-01T00:00:00.000Z', '2444-01-01T00:00:00.000Z', '2463-01-01T00:00:00.000Z', '2475-01-01T00:00:00.000Z', '2489-01-01T00:00:00.000Z', '2500-01-01T00:00:00.000Z', '2518-01-01T00:00:00.000Z', '2535-01-01T00:00:00.000Z', '2545-01-01T00:00:00.000Z']You can see that we have the times displayed here in the same format as we have used before (tweaked time values, not the geological age in millions of years). Now that we have a list with all available ages, it is possible to iterate through this list and analyze the map associated with every time step.
In the code snippet below, we create a for loop that will use each value in the times list, extract the geological age value, create a WCS request and get the data, calculate its land and ocean area, and store them in a dictionary.
results = {}
for time in times:
geological_age = int(time[:4]) - 2000
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}"
)
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,
False
)
results[geological_age] = {
'land_area': float(land_area),
'ocean_area': float(ocean_area)
}
del raster, data
print(results){0: {'land_area': 158193714972089.66, 'ocean_area': 351871900996932.5}, 6: {'land_area': 155229003015959.22, 'ocean_area': 354836612953062.9}, 11: {'land_area': 152970779439717.3, 'ocean_area': 357094836529304.8}, 15: {'land_area': 156210409047566.06, 'ocean_area': 353855206921456.06}, 20: {'land_area': 154818326029319.66, 'ocean_area': 355247289939702.5}, 33: {'land_area': 152344567014130.56, 'ocean_area': 357721048954891.56}, 40: {'land_area': 148632678865657.3, 'ocean_area': 361432937103364.8}, 48: {'land_area': 146346163763802.22, 'ocean_area': 363719452205219.9}, 56: {'land_area': 140238893123393.42, 'ocean_area': 369826722845628.7}, 68: {'land_area': 135971671318449.17, 'ocean_area': 374093944650572.94}, 84: {'land_area': 122434026357452.64, 'ocean_area': 387631589611569.5}, 94: {'land_area': 129162011070413.84, 'ocean_area': 380903604898608.3}, 100: {'land_area': 135611979059946.42, 'ocean_area': 374453636909075.7}, 113: {'land_area': 94771512331389.27, 'ocean_area': 415294103637632.9}, 120: {'land_area': 135358854880255.23, 'ocean_area': 374706761088766.9}, 133: {'land_area': 143445632583091.12, 'ocean_area': 366619983385931.0}, 140: {'land_area': 159986877851679.3, 'ocean_area': 350078738117342.8}, 154: {'land_area': 159596894666459.97, 'ocean_area': 350468721302562.1}, 165: {'land_area': 168327179614427.12, 'ocean_area': 341738436354595.0}, 180: {'land_area': 164209413041962.47, 'ocean_area': 345856202927059.7}, 200: {'land_area': 145695358704643.56, 'ocean_area': 364370257264378.56}, 210: {'land_area': 142902195363185.38, 'ocean_area': 367163420605836.75}, 220: {'land_area': 143147821788738.8, 'ocean_area': 366917794180283.3}, 230: {'land_area': 143492018688690.16, 'ocean_area': 366573597280331.94}, 240: {'land_area': 150918594147181.4, 'ocean_area': 359147021821840.7}, 250: {'land_area': 153456433967734.66, 'ocean_area': 356609182001287.5}, 270: {'land_area': 157502821920807.72, 'ocean_area': 352562794048214.44}, 290: {'land_area': 159834823397765.28, 'ocean_area': 350230792571256.8}, 300: {'land_area': 165935895894753.53, 'ocean_area': 344129720074268.56}, 315: {'land_area': 158615988485129.0, 'ocean_area': 351449627483893.1}, 331: {'land_area': 157735752149354.5, 'ocean_area': 352329863819667.6}, 350: {'land_area': 149498419543432.12, 'ocean_area': 360567196425590.0}, 370: {'land_area': 145220500942584.66, 'ocean_area': 364845115026437.5}, 383: {'land_area': 144787030783366.22, 'ocean_area': 365278585185655.9}, 393: {'land_area': 144431237357015.12, 'ocean_area': 365634378612007.0}, 408: {'land_area': 147071746424228.25, 'ocean_area': 362993869544793.9}, 420: {'land_area': 149528310589928.03, 'ocean_area': 360537305379094.06}, 444: {'land_area': 139432934538610.47, 'ocean_area': 370632681430411.7}, 463: {'land_area': 139693856382604.95, 'ocean_area': 370371759586417.2}, 475: {'land_area': 139471622949961.38, 'ocean_area': 370593993019060.75}, 489: {'land_area': 141579691503337.5, 'ocean_area': 368485924465684.6}, 500: {'land_area': 140892897224317.53, 'ocean_area': 369172718744704.6}, 518: {'land_area': 142827817642138.7, 'ocean_area': 367237798326883.44}, 535: {'land_area': 151802129494776.53, 'ocean_area': 358263486474245.6}, 545: {'land_area': 137047848962357.45, 'ocean_area': 373017767006664.7}}Finally, we can plot the evolution of these two variables through time. In the example below, we create two separate plots, one for land, another for oceans with a reversed x-axis value (0Ma = present-day is on the right-side).
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=(8, 3.5))
ax1.plot(ages, land_areas, color='darkgreen', marker='o', markersize=4)
ax1.set_xlabel('Age (Ma)')
ax1.set_ylabel('Land area (m²)')
ax1.invert_xaxis()
ax1.set_box_aspect(1)
ax2.plot(ages, ocean_areas, color='darkblue', marker='o', markersize=4)
ax2.set_xlabel('Age (Ma)')
ax2.set_ylabel('Ocean area (m²)')
ax2.invert_xaxis()
ax2.set_box_aspect(1)
plt.tight_layout()
plt.show()
→ What do you observe ?
→ You have seen what we can do with simple masks representing land and ocean. Do you have any ideas about other statistics we could easily plot ? Try it !
transform attribute. What does this attribute describe ?Web Map Service (WMS) — An OGC standard protocol for requesting georeferenced map images (e.g., PNG, JPEG) from remote servers.
Web Coverage Service (WCS) — An OGC standard for accessing raw raster datasets (“coverages”) such as elevation models, climate grids, and satellite imagery, preserving full numeric values.
matplotlib A comprehensive plotting library for creating static, animated, and interactive visualizations in Python.
numpy The fundamental package for efficient numerical computing in Python, providing fast arrays and mathematical operations.
rasterio A Python library for reading, writing, and manipulating geospatial raster data using GDAL under the hood.
requests A simple and elegant HTTP library for sending web requests and interacting with APIs.
xml.etree.ElementTree A Python standard library module for parsing and navigating XML documents as a tree structure, where each element’s tag, attributes, and text are accessible as nodes.