Skip to content

Cartography as code: Cartopy

This story calls back a long way to my mapping roots. I started out making scientific maps using the venerable Generic Mapping Tools (GMT) to plot bathymetry, ship tracks, sampling stations. Almost two decades later, I’m still making maps with code! However, this time it’s in Python – with some helpers.

Here’s the end product – a timelapse video of sea ice moving across the Arctic Ocean and the Barents Sea, showing daily states of the sea ice edge, sea ice concentration and it’s velocity across the surface of the ocean. It was assembled in Python from various data sources, with helpers along the way to do some cartographic parts and some final video assembly. The moving grey line shows the sea ice edge. Ice concentration is visualised with a blue to white gradient (white is 100% sea ice cover), and the grey squares show ice velocity (darker grey is faster).

A caveat to this post: work is underway on the code behind this movie – I reserve the right to make it better and update at any time!

The brief

The brief for this map was fairly simple. While I had a lot of Antarctic experience, when I arrived in Norway in 2020 my knowledge of Barents Sea ice dynamics was slim. So I needed a matrix-style fast upload. I had a lot of static graphs at my disposal, however I needed more – so I designed a simple map to animate daily sea ice edges. This was the key! So I kept going.

What I discovered I needed was a way to look at historical ice motion and dynamics over the Barents Sea. And ice inflow from the Arctic Ocean. And also show how ice was present or absent in different seasons or years. And also show the fraction of ice covered by ocean. Once I started and realised it was useful, it became a tool to deliver at workshops and conferences. So the brief became:

  • Display daily sea ice concentration
  • Display daily sea ice edge
  • Visualise ice motion
  • Add geographic context
  • Within institutional design guidelines, display project and institutional branding.
  • Assemble a timelapse movie for a few years of daily data.

Here’s how it was met.

Making a basemap

Creating the ocean floor basemap is a story almost of its own about terrain rendering methods. In this region the primary issue is that we want to display regions of extreme relied – the four thousand metre deep Arctic Basin and thousands of metres tall Gakkel Ridge, on the same map as some shallow, very low relief regions in the Barents Sea – just a a few hundred metres deep in total.

The next issue is that no one dataset covers the whole region. I used IBCAOv4 and GEBCO bathymetry, merging a subset of the North Atlantic with the Barents, creating a kind of keyhole shaped map to ensure that all the regions I wanted to show in the end product had some map coverage.

The rendering stack is fairly straightforward. Using an inverted layer order (1 is the bottom), the stack is:

  1. IBCAOv4 + GEBCO gridded terrain, tinted with a blue (deepest) to white (sea level) color ramp, blend mode normal
  2. A unidirectional hillshaded copy of the same terrain, vertical exaggeration (Z-factor) set to 7, blend mode multiply
  3. A multidirectional hillshaded grid, vertical exaggeration set to 2, blend mode multiply
  4. Slope derived from the same grid, black (steepest) -> white (shallowest) colour ramp, blend mode overlay
  5. Hillshaded slope, blend mode burn
  6. Optionally to make everything really contrasty, another hillshaded slope layer with blend mode overlay.

I settled on using only layers 1-5, then rendered out a complete shaded geotiff to ingest into Python. It took a few rendering layers to get just right – mainly getting some features popping up in the Barents Sea. Fun fact: my 2023 Geohipster calendar map is a subset of this baemap!

Choosing data sources

In sea ice research, choosing the large scale datasets used to derive insights is a political act. In one corner is the NASA team. In a second corner is a strong research group at the University of Bremen. And a third contender is the Copernicus Marine program (OSISAF). There are endless publications about whos data are better or not. My choices were driven by two things: First, I wanted the data to look like what we see in the real world. Second, I wanted a specific visual approach. In the end all three data providers played a role!

NASA’s MASIE sea ice edges were chosen to represent sea ice edge because they are drawn from multiple sensor data and were available continuously for the time period.

Sea ice concentration was visualised using the University of Bremen’s passive microwave dataset. Again, it was continuous for the timeframe. It was also ‘just high enough resolution’ at 6.25km a pixel. This means it looked more like ‘reality’ than an over-smooth larger scale dataset. Even at a whole of Arctic scale, showing a 25 or 30km per pixel dataset leads to a false impression of sea ice dynamic states!

…and now I’ll eat those words because I used a 30km per pixel sea ice motion product from OSISAF. My rationa.. actually let’s be accurate: excuse is that ice movement isn’t as critical to show at fine scales.

Data from all 3 sources was mostly continuous for the whole time series. Because gaps appear in different datasets at different times, an early task in processing is to ensure that a daily plot doesn’t get out of sync – for example if there are four days of ice concentration missing, we need to ensure that those days are handled properly. You’ll see that gaps exist in the ice motion data relatively often early in the time series.

You’ll see in the video some motion days are missing – on these days there was a dataset there, it just didn’t contain valid data.

Making the movie map

Because I needed to automate making hundreds of maps, I chose Python as a cartographic tool. The map is assembled in a Jupyer notebook which also takes care of data filtering and preprocessing. Here is the notebook: https://gitlab.com/adamsteer/aen/-/blob/main/jupyter-notebooks/ice-conc-drift-nearsideperspective.ipynb

I used Cartopy as the map renderer, because it integrates well with the familiar Matplotlib interface – adding a projection handling capability and very conveniently, a Natural Earth data import handler for coastlines!

Setting up the processing environment

I use the conda package manager on linux for Python work. And tend to do things using Jupyter notebooks. After installing miniconda, the next steps are:

  1. create a virtual environment for generalised Jupyter notebook Cartopy work
  2. start the environment
  3. start Jupyter

Just under the block where all the library imports are done, you’ll see a cell starting with a bang !. This is a way of Jupyter interacting with the system to do things – so, if you get complaints about missing Python libraries, you can grab them and restart:

! conda install [library] -y

My environment looks like:

conda create -n jupyter jupyter cartopy rasterio netcdf4 h5netcdf xarray os pandas -c conda-forge

…and the notebook is run like:

(base)#> conda activate jupyter
(jupyter)#> jupyter notebook

Now, time to make the map!

Styling decisions

Deciding how to style the sea ice was a process. We quite often see sea ice rendered in garish colours to differentiate it from the ocean and other things. I went the other way – using a very similar colour ramp to the basemap. My design choices were driven by the fact that, well, sea ice is bright, and places with less sea ice are in fact blue (or ocean colour – not necessarily just blue). After some experimentation and asking colleagues for review, I settled on the white (lots of ice) to blue (no ice) scale. It is easy to look at and nobody mistook icy areas for ocean floor. Being able to still see the ocean floor below the ice is appealing to me – it gives a good indication about any topographic associations with sea ice motion. Really, it’s the topography pushing water which pushes ice, although in some cases we get these huge open areas leeward of islands – where accelerating wind pushes ice offshore faster than it gets made. This isn’t the place for a full discussion of what drives ice dynamics – for a generalised explanation it’s enough to understand that ice motion is the integral of wind and water influences each of which varies. Topography can influence both water and wind, thereby influencing ice motion.

The MASIE ice edge is a grey line – again chosen because while less distinct, it is easy (for me) to look at. It could be more distinct, however it probably not the most important visual element. From my view, ice concentration is the primary thing to look at because it tells us more about ice states than just an edge of ice cover.

Perhaps the most complex style decision was ice motion. I experimented with drawing arrows and other vector means. While informative, it was just…ugly. In order to be visually satisfied I had to reduce the data resolution too far. The perils of science design! Balancing a need for accuracy with a need to be able to look at it. I was really just tweaking knobs when the idea of using shading to show motion popped up. After quite a few renderings of short videos I thought that shading grid cells to show motion worked in terms of showing drift motion patterns, and also added an interesting visual element to the map.

Here is a single frame example:

Animation

Originally this map was animated in Python. This was risky – because it took a really long time. It also meant that changes to frame rates and other tweaks meant starting over. I settled on having Python produce both a video and daily frames, then assembling a timelapse using ffmpeg. There are many many instructions for how to do this on the internet!

Whats wrong with this map?

In the movie you’ll see that there are ‘ice motion’ data very far south of actual ice – a little ring around the Arctic circle. This is a cautionary tale about over-reliance on numbers from remotely sensed data without mapping them out – or heading into the field to validate. It is not a criticism of a research group or approach, it just is. Even using surveying GNSS and drones and digging snow pits by hand, we get data artefacts that must be checked before use.

The labels could also be better – I really should add labels for Russia, Nova Zemlya, the Svalbard archipelago. And the Barents Sea, and the Arctic Basin. One day perhaps.

I think the ice motion data is sometimes just too jumpy, I did try ffmpeg’s interpolation to fix that, however it made a total mess of the datestamps.

What can we learn from this movie?

It’s just another sea ice motion animation yeah? Well, yes. I was extremely interested in what happened in the Barents Sea – where ice came from, how much, what did it look like outside of line charts showing areas. What I got out of this was an understanding of how dynamic the ice situation is in the Northern Barents. It is important to understand whether ice grows there or if it is imported – and it looks a lot like heavy ice years in the Barents sea are associated with a lot of motion and some cycles of northerly winds pushing ice south from the Arctic. In the field it was sometimes pretty obvious to spot Barents ice and Arctic Ocean ice, however not always…

I also used this movie as a prop in talks to the wider project community about sea ice dynamics. It helped a lot to show how, if we travelled into the field in November and March, the ice edge might have retreated but also it wasn’t linear – there are pulses and eddies and many other things happening. So if we write a paper that says ‘the ice edge was at X in Spring and Y in Summer’, we can expand now to say ‘the ice idge was at X in spring and Y in summer, with [ retreat | advance ] impacting light availability in the ocean between sampling stations’. So the inference that sea ice comes and goes more or less linearly can be put to bed properly. Yes, it kinda does – it also doesn’t in ways that impact the system.

Sub learning: Why is MASIE ice edge data not exactly on University of Bremen ice edges?

You should notice that ice edges (grey lines) are not always aligned with the edge of ice concentrations. This is because of one of two things:

1. each uses a different suite of sensors and data resolutions to assess whether sea ice is present on the ocean or not
2. I’ve made a mistake in my code and offset days accidentally

Even if item 2 is not true, item 1 is always true and highlights a reason to use multiple sensors and products whenever making assertions about earth systems. Each sensor has an implicit bias, enhanced by the code used to generate a product. It is extremely important to note that each product is made with accuracy as the first and foremost concern. No researcher wants to create misleading data, everyone does the best they can with the tools they have. If you want to split science hairs, my suggested approach is to look at the story rather than the words, to use a metaphor. Do all the datasets tell a common story? Great! We’re good. It doesn’t matter that one line is drawn slightly different to the other – as much as careers are made on these things.

It is also extremely difficult to validate the datasets used here – when we go into the field we observe scales of metres to maybe a kilometre or two. These datasets measure pixels with multiple kilometres per side. Scale matching between what we measure on ground or view from a ship and satellites ain’t easy. I’ve got at least a decade of working toward that goal under my belt.

Sea ice is complex. And never stays where you left it…

Summary

This is a movie showing sea ice dynamics in the Arctic mapped from public open data using Cartopy in a Jupyer notebook to assemble it all. It also walks through a technique for displaying mixed-relief topography in QGIS. The end result was both a tool for my own understanding, and a valuable prop for discussion around sea ice dynamics and framing research outputs with accurate concepts about sea ice states in mind.

I think we understate this type of work in research. It doesn’t have a paper in mind, it was never aimed to be published. It’s value lies outside of the standard value system of research. Yet it is important and valuable work, and the kind of thing we need to do often if we want to create publications which claim expertise on a topic. A criticism I still hold for much research is that we rely on static graphs and derived datastets 3 or 4 or more steps removed from the on-ground reality. If our work is determined by ‘is it going in a paper/how many papers can we get from this?’, we all lose.

A solution I propose – and have staked my career on – is investing in the kind of work you see here to aid in system understanding. And also spending hours watching ice from the lowest deck of the ship instead of writing papers in a cabin. Researchers need the freedom to do these exploratory works, take an artists mindset, become absorbed. We celebrate the heck out of people who did it in the past…

Take a look at the notebook: https://gitlab.com/adamsteer/aen/-/blob/main/jupyter-notebooks/ice-conc-drift-nearsideperspective.ipynb for how everything was assembled in Cartopy – it walks through quite a bit about projections and views and Z ordering of graphical elements. I’m totally happy to take issues (don’t hold your breath awaiting fixes) and pull requests.

Above all, I hope all this is useful to you – or at least an interesting read.

The sales pitch

Spatialised is a fully independent consulting business. Everything you see here is open for you to use and reuse, without ads. WordPress sets a few cookies for statistics aggregation, we use those to see how many visitors we got and how popular things are.

If you find the content here useful to your business or research or billion dollar startup idea, you can support production of ideas and open source geo-recipes via Paypal, hire me to do stuff; or hire me to talk about stuff.