Motivation

A team in the Informatics Lab is working on applying machine learning to the thorny question of how to predict weather in a region of the world where seasonal weather forecasts have low skill. The region in question is North East China, where 5% of global maize is grown. Weather in the region is heavily influenced by monsoons and weather year-to-year can be highly variable, impacting the success of maize growth and harvests.

The project depends upon separate activities written in both Python and R, reflecting its diverse requirements and interdisciplinary nature. The purpose of this blogpost is to share our experiences of making analysis codes openly available and reproducible by the public, a key requirement underpinning open science. It will discuss how to configure a Binder to handle both R and Python code, using literate programming to encapsulate each analysis with a narrative justifying its purpose in the overall research.

Prior art

Many will be familiar with Jupyter Notebooks, and by extension Binder. For the uninitiated, Jupyter Notebooks are documents containing text, images and executable code which explains some sort of analysis with a narrative and working code. It is a front end to a live backend which executes the code cells in the selected language (Python, R, Julia, etc.).

Binder is an automatic service for launching a live, executable Jupyter Notebook (and other frontends) on the cloud, with the infrastructure gracefully handled by Docker. The generous people at Binder provide free, limited compute resources for the purposes of sharing research, demonstrators, educational material and even documentation. This cloud computing is provided for free by supporting organisations such as Google Open Source and the Alan Turing Institute. (The instance used in this example and blogpost is run on the “default” Google Cloud infrastructure). The whole project is ran by volunteers, and financially supported by donations and grants. The underlying technology can even be deployed on your own (organisation’s) infrastructure. Usage guidelines can be found here.

From a user perspective, it works by simply pointing the service at an online git repository (e.g. hosted on GitHub) and Binder handles the rest. It has been around for a few years and has been continually developed to now include other user interfaces such as Jupyter Lab (an IDE like implementation of Jupyter Notebook), RStudio and RShiny.

RStudio will be familiar to those who write R code as the go-to IDE for developing R code. Traditionally used as a desktop application, RStudio now supports being run in a browser with RStudio Server and is available as a service through RStudio Cloud.

For this project I followed many guides on how to deploy a Binder for both R and Python code, since (maybe unsurprisingly) there is more than one way to configure the infrastructure of a Binder. I started with Zero-to-Binder, the canonical guide hosted by The Alan Turing Institute as part of The Turing Way. While this is a very useful guide, and setup for Binder with Python “Just Works™️”, I found that its recommendations around setting up a Binder for R were not as functional as implied. Following the method for R outlined in the main guide leads to a working Binder for (just) R. However, it didn’t support a Python Binder as well and I was not satisfied with how the R package installation was handled. So I added renv to the postBuild script to handle package specification more explicitly (see Working Configuration below for more details) and searched for some guidance on Binderising R and Python together.

As an aside, the guide recommends using the holepunch package, an R package which automatically configures a Binder instance for your R code. Unfortunately, it has not been maintained for over two years and I found it produced a completely blank Dockerfile that did not produce a working Binder instance, so I abandoned that line of action. The guide also recommends using Rocker, a repository of R-specific Docker images on DockerHub. While this works, this level of manual configuration is unnecessary with the now excellent automagic of repo2docker, so I would not recommend using this method for setting up a Binder either (although very useful for setting up a remote RStudio/RShiny instance!).

I ended up taking inspiration from the binder-examples/r_with_python and binder-examples/r repositories, which led to the final configuration described below.

Working configuration

The configuration that I found to work well for RStudio, Jupyter Lab and RShiny is demonstrated in this GitHub repository:
https://github.com/informatics-lab/binder_rstudio_jupyterlab_example

The infrastructural specifications are found in the .binder/ folder and consist of:

.binder/
├─ environment.yml   # Define *additional* conda packages to install
├─ install.R         # Define R packages to install (not used)
├─ postBuild         # Where renv::restore() is executed
├─ renv.lock         # Manifest of R packages to install using renv
├─ runtime.txt       # R version (from the MRAN archive)

• environment.yml is used by conda to install additional packages for Python, including from PyPI using pip. Binder automatically includes a bunch of base Python packages such at jupyter, python, IPython, etc, so you do not need to define these yourself. It is possible to install R packages here too but not recommended as the conda channels for R are not well maintained.
• install.R is the “standard” way of installing R packages for a Binder using install.packages() or devtools::install_version(), but I chose to use renv instead since it is more reproducible and allows you to explicitly specify package versions and sources in a single manifest file, renv.lock.
• postBuild is run after the Docker image has been built, and this is where renv::restore() is executed. There is discussion about integrating renv into the Binder build process but that hasn’t been implemented so for now this is the easiest solution. It adds a little overhead to the Binder spin-up time, since the renv environment is not cached with the Docker image, but it is fast enough in my experience.
• renv.lock defines the list of R packages to install using renv::restore().
• runtime.txt is used to define the version on R (from MRAN) to install. It uses the format r-<version>-<year>-<month>-<day> as the only line of the file (I tried adding comments but it very broke the system)

This relatively minimal configuration setup allows you to deploy a Binder with Juptyer Lab, Jupyter Notebook, RStudio and RShiny as the interfaces, all without having to define or edit a Dockerfile. This is all down to the magic of repo2docker and binderhub — kudos to the JupyterHub team!*

*The ease of setting up such a system might raise the question of how secure Binder is. The Binder team have take precautions, as discussed in their docs.

Examples

There are four demonstration analyses included in the binder_rstudio_jupyterlab_example repository. All four plot example data for maize yields in three Chinese provinces, acquired from the National Bureau of Statistics of China, and ERA5 reanalysis weather data for those same provinces.

Jupyter(Lab) with Python

The orange JUPYTERLAB badge in the link above will open a Jupyter Lab session with the Jupyter Notebook analysis-4.ipynb open. This notebook plots data in the data/ directory as scatter plots with linear regression lines fitted using Python. You can execute the cells in the notebook which will be run using the live Python kernel in the background. You can add code, edit the notebook, create new notebooks and files, and add your own data to analyse (using the file browser on the left hand pane).

Jupyter(Notebook) with R

The orange JUPYTERNOTEBOOK badge in the link above will open the Jupyter Notebook analysis-2.ipynb in a standalone Jupyter Notebook session. You can run the cells like a regular Jupyter Notebook, which plots data in the data/ directory as scatter plots with linear regression lines fitted, using an R kernel. You can add code and edit the notebook but you cannot easily add data to analyse. The JupyterLab example above has a better user interface for adding and navigating files. This usecase is not specific to an R-based Jupyter Notebook, simply a demonstration of another interface.

RStudio

The blue RSTUDIO badge in the link above will launch an RStudio session, from which you can navigate the file browser to open the R script analysis-1.R and RMarkdown notebook analysis-3.Rmd. These can be run as if you were using RStudio Desktop. You can edit these scripts, create new files, and upload your own data for analysis.

RMarkdown notebook in RShiny

The blue RMARKDOWN badge in the link above will open the RMarkdown notebook analysis-3.Rmd in an interactive RShiny instance. It simply plots data in the data/ directory as scatter plots with linear regression lines fitted, and embeds a couple of RShiny widgets. Cells cannot be executed or edited like a Jupyter notebook, but the widgets can be interacted with like in an RShiny dashboard.

This configuration took a lot of digging to uncover an implementation for. The functionality is well established within in R, but its implementation within the Binder project was unclear. These two pull requests (#799 & #891) in the repo2docker GitHub repo eventually led to a solution. Please note that in order for this setup to work the .Rmd file must be located in the base directory of the repo, subfolders are not supported for this functionality. The name of the file does not matter.

RShiny dashboards are also supported in MyBinder using a server.R and ui.R configuration. This setup is not demonstrated in this repo.

Future work

As a result of doing this project there are a couple of things that could improve the system:

• Saving the final iteration of the Binder Docker image out to a container registry like DockerHub, where fully built images are saved and then downloaded and run for quick spin-up times (since the building stage has already been done). There does not seem to be an easy way to do this with repo2docker or binderhub, but would be useful for an organisation like the Met Office to provide pre-compiled Docker images for the finalised software environments used to recreate analyses for specific papers/projects.
• Include renv in the R Binder stack. This is discussed in the following GitHub issue on the repo2binder repository.
• Support for RMarkdown notebooks could be better. At present the implementation seen in the example above requires that the .Rmd file be located in the base directory of the repository. Support for locating the files in a subdirectory would allow more flexibility for repository configurations.
• Using flexdashboard to arrange RShiny dashboards. It looks like a cool technology that builds on the ability of RMarkdown to be used with RShiny.

Conclusion

Having been aware of and used Binder for over four years now, it’s exciting to see it go from strength to strength. It’s an invaluable resource for creating demonstrations and reduces the barrier to reproducible research. While the configurations demonstrated in this blogpost took a bit of digging, trial and error, it was a surprisingly painless (and fun) week where very little broke completely or just didn’t work. Kudos to the community and Jupyter devs for working so in the open that there are plenty of breadcrumbs to follow for an answer.

Given the continued ease of setting up a Binder instance, I encourage all scientists, project maintainers and research software engineers to have a go with Binder — it might be just the (free) resource you’ve been looking for.

R and Python got you in a Bind? was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

A New Order

November 2020 saw the announcement of a significant change of direction for Apple: for the first time in nearly 15 years, the Apple Mac computer would no longer use Intel processors. Intel processors would be replaced by the M1 System on Chip (SoC), developed in-house by Apple and an enhancement on the A-series SoCs that have powered iPhones and iPads for a number of years.

Why does this matter?

Intel’s x86 microarchitecture is very prevalent across desktop computers, laptops, servers, the Cloud, and HPC (although IBM’s Power microarchitecture is a notable exception in the HPC sphere). ARM’s microarchitecture, meanwhile, is completely ubiquitous in mobile devices. If you own a modern smartphone, it’s almost certain to be running an ARM-architected chip, just manufactured under license by Qualcomm, Samsung, Huawei, MediaTek, Apple (of course), or another chip manufacturer.

The move by Apple to using the M1 chip for desktop computer applications is much more significant than turning MacBooks into glorified mobile phones, however. The change in microarchitecture from x86 to ARM introduces new uncertainty to the usability of desktop applications. Putting aside the question of the different application installers used by Windows, macOS and Linux, the likelihood was that if you compiled a binary (that is, a program compiled to an executable object) for x86, it would be likely to function on all of these desktop operating systems without trouble.

The M1 SoC’s ARM microarchitecture upends this comfortable certainty. A different microarchitecture means a different instruction set, which fundamentally means that a binary compiled for x86 will not run on ARM. It’s similar to going into a café in a country where you don’t speak the local language — it’s going to be hard to ask to buy a local delicacy if you and the barista don’t understand each others’ spoken language.

The Pertinent Questions

To get around this, you now require all the programs that you want to use on an M1 to be compiled for ARM microarchitecture as well as x86. The question is, have they been?

This question recently came to roost in the Informatics Lab as we received our first Macs that use the M1 SoC. The Lab (and indeed the wider Met Office) run a lot of Python-based Data Science and Machine Learning analysis workflows on large n-dimensional gridded weather and climate datasets.

The pertinent question for the Lab, therefore, is can we run these workflows on an M1 Mac? And if we can, how does the experience compare to running the same workflows on an Intel Mac? We’ll spend the rest of this blog post exploring the answers to these questions.

Data Science

I’ll handle Data Science as distinct from Machine Learning here and define it as all classical analysis workflows that don’t involve an element of machine-driven learning (unsurprisingly we’ll look at that in the next section on Machine Learning).

The Lab’s data science workflows typically use common scientific Python packages to load, analyse, visualise and save nD datasets of numerical data describing weather and climate phenomena. These are often provided as CSV or NetCDF files. To do this we use Python libraries such as Jupyter, NumPy, matplotlib, dask, pandas, Iris and cartopy, and Xarray, and their surrounding ecosystems of dependency packages. The simple question I wanted to ask of the M1 Mac in this case is “can I install and use these Python libraries on an M1 Mac?”

To test this I installed Python and conda via miniforge and went to work setting up environments that contained these different packages. Quite simply, I found that all of these packages installed and functioned without issue on an M1 Mac, so this was an easy and very gratifying conclusion to reach.

There is no doubt that this state of affairs is on account of the amazing efforts put in by the maintainers — both of the libraries themselves and the package builds of the libraries at conda-forge. Their efforts have ensured that these libraries are all available with either microarchitecture-independent (noarch) package builds, or with builds for multiple different microarchitectures, including x86 and ARM. So, when you install these libraries using conda, it will automatically select the package build that’s best suited to your system, ensuring consistent functionality across all sorts of different systems.

Machine Learning

For Machine Learning (ML) workflows I wanted to take things beyond the first question of “can I install the libraries?” and explore whether we can also find a computational advantage for using the M1 Intel. There is good reason for hoping that we might: the M1 Pro SoC in our MacBooks Pros contains an 8 core CPU as well as a 14 core GPU and a 16 core Neural Engine. These last two items are both accelerators that can be used to speed up the training of ML models — a GPU handles matrix operations at high speed typically to drive graphical applications, but can just as well be applied to matrix operations in ML applications; and the Apple Neural Engine is a type of Neural Processing Unit (NPU) that’s specifically designed for the sort of operations found in Neural Networks.

However, just because an accelerator is present doesn’t necessarily mean that you can just use it to accelerate training your ML model. You need to have a mechanism to transfer the training to the accelerator otherwise it will just run on the CPU and you won’t benefit from the presence of the accelerator. This fact turned out to particularly apply to the Apple Neural Engine, which seems to only be accessible when using Core ML, Apple’s Swift library for running ML models. But as this is written in Swift it’s not a lot of direct use for Python ML libraries, so I didn’t explore this further.

The GPU on the M1 SoC, however, is accessible via Metal, Apple’s graphical accelerator API. Apple have also written a TensorFlow plugin for Metal that enables accelerated training of ML models written in TensorFlow using Mac GPUs. It made sense, therefore, to test accelerating ML workflows using ML models built with TensorFlow (PyTorch and XGBoost, for example, only support CUDA-based acceleration on NVIDIA accelerators).

The Experiment — TensorFlow

Accelerators work best when data volumes are large and the computations being run are complex. This is because there is a computational overhead associated with moving the model and in-memory data from the CPU to the accelerator. If the model and/or data are small, this overhead may be significantly larger than any speed-up gained from running on the accelerator, and the net result is that performance actually decreases for using the accelerator.

So, to best test for performance improvements for using the GPU in the M1 SoC, we want a complex ML model to be run on large data volumes. Or perhaps better still, we want a series of models (ranging from simple to complex) but all running on the same data, from which we can look for noticeable performance improvements when accelerating the more complex models compared to the simpler models.

Pleasingly, just such an example can be found in the Overfitting and Underfitting tutorial in the TensorFlow docs. This demonstrates an important principle in ML model architecture (a complex model doesn’t necessarily produce better results) by running a series of models of increasing complexity against the same large-ish dataset. Or, to put it more simply, this is exactly what we were after. More details about the architectures of the different models and the dataset used can be found at the link above.

The setup of the experiment was as follows:

• Use six of the models defined in the tutorial linked above (the models named ‘tiny’, ‘medium’, ‘large’, ‘L2’, ‘dropout’ and ‘combined’). Record the wall time taken to train the model and the total number of training epochs required, and find the result value of number of seconds per epoch for the training run.
• Each model was run three times on each hardware platform then the mean result value across the three runs was calculated.
• As Metal allows us to accelerate ML models using any recent Mac with a GPU, we can compare the M1 SoC to the previous generation of Intel+AMD Macs.

Therefore, the following hardware platforms were used:

• A 14” MacBook Pro 2021 with a M1 Pro SoC with an 8 core CPU, a 14 core GPU and 16GB integrated memory, with and without Metal acceleration.
• A 16” MacBook Pro 2019 with an 8 core Intel i9 9880H, an AMD Radeon Pro 5500M 4GB and 32GB system memory, with and without Metal acceleration.

Thus we had two devices each running with and without Metal-enabled acceleration, providing four platforms in total, on which to run six experiments three times each. To accurately compare different models that ran to a different number of epochs, the data values we chose to record was the time in seconds taken to calculate one training epoch, calculated by dividing the total number of wall-time seconds that the model trained for by the number of epochs taken to train.

Results — TensorFlow

The graph below shows the average runtime (s) per epoch for each of the six models on each of the four hardware platforms.

The graph demonstrates a number of interesting results:

• The M1 SoC outperforms the Intel processor both with and without Metal acceleration enabled.
• Metal acceleration makes small and medium sized ML models a lot slower to train than without Metal acceleration.
• However, as the models get larger, the difference between Metal acceleration being enabled and not being enabled shrinks for the M1 Mac, and in fact the Intel Mac trains faster for the two largest models with Metal acceleration enabled.
• The performance difference between M1 and Intel gets more pronounced as the ML models get more complex.

One other result that is not demonstrated in the graph above is the thermal impact of running these experiments, particularly with the most complex models and with Metal acceleration enabled. The M1 Mac became warm to the touch, but at no point did its fans switch on. The Intel Mac, meanwhile, became too hot to touch and its fans were running at maximum — to the extent that some runs showed obvious evidence of thermal throttling.

The Experiment — PyTorch

PyTorch has just (at time of writing) released support for Metal devices in v1.12. Unlike TensorFlow’s Metal plugin which is provided by Apple, the PyTorch plugin has been provided by the PyTorch community. This means we can run a similar experiment in PyTorch to the experiment we ran in TensorFlow.

Again we’ll use an example from the documentation to provide a reasonably-sized ML model with which we can test the performance of Metal acceleration. For PyTorch we’ll use this transfer learning tutorial and use both the fine-tuning and feature extraction applications. We’ll use the same hardware platforms for PyTorch as were used in the TensorFlow experiment.

Results — PyTorch

The bar chart below shows the results of running the two applications from the transfer learning tutorial on the Intel and M1 Mac, with and without Metal acceleration enabled. Here we simply recorded the total wall time for each application on each platform and treated these as our results.

Again, we see broadly the same pattern as with the TensorFlow experiment — typically the M1 is faster than the Intel Mac, and there is typically limited benefit for using the GPU (if any speed-up at all). The notable exception is with the fine-tuning experiment on the M1 Mac, where a significant speed-up was recorded with the GPU enabled.

Note that the results for this experiment were recorded a little differently than to the TensorFlow experiment, with total time elapsed for each application recorded rather than time per epoch. Each experiment was also run only once, meaning spurious results will not be averaged out over multiple runs.

Timing each run allowed us to capture not only wall time (the time between start and end of execution) but also the user process time (the amount of time the CPU spent executing the user process). CPU time for the two applications, with and without GPU acceleration, are shown below:

It’s significant that CPU time decreases significantly when the GPU is enabled. This is a direct demonstration that the processing is being offloaded from the CPU to the GPU, so the time the CPU spends running the processing is significantly reduced, even if the overall wall time of the processing (from the previous chart) does not change significantly.

Conclusions

The results from both these experiments allow us to draw these interesting conclusions:

• There is an evident improvement in training more complex ML models on an Intel Mac with Metal acceleration enabled.
• There is the suggestion of a similar performance improvement on the M1 Mac if the ML model was even more complex.
• The difference in thermal performance between M1 and Intel is quite striking: the M1 is not only faster but also far more thermally efficient than Intel. Indeed this conclusion reflects one of the main reasons that Apple pursued an ARM-based SoC for Mac computers.

In Summary

This brief exploration of using an M1 SoC for Data Science and Machine Learning workflows has found a number of conclusions that are significant to the applicability of the M1 to these fields:

• We encountered no difficulty with installing and using common Scientific Python libraries for Data Science workflows on an M1 Mac.
• The M1 SoC is faster than Intel in the ML model training experiment that we performed. And the M1 seemed to be faster for the Data Science workflows that I used to test the functionality of scientific Python libraries on the M1 Mac.
• We found that Apple’s Metal API can be effectively used to accelerate training ML models both on M1 and Intel Macs, with a performance gain evident with more complex models on the Intel Mac.
• The M1 SoC is much more thermally efficient than Intel.

Next Steps

Here are some ideas of where to go to extend this investigation further:

• Data Science workflows felt like they were faster on the M1. We could time these workflows and see if we can prove (or disprove!) this with some numerical results.
• The ML experiment showed that the M1 was always faster without Metal acceleration. It would be interesting to explore how big an ML model we would need for the M1 with Metal acceleration to train faster than without.
• I still want to try and make use of the Neural Engine to accelerate ML workflows. It would be a lot of effort, but it would be interesting to see if we could reproduce the ML model in Swift for Core ML, or even call Core ML from Python somehow and see if we could make use of the Neural Engine.
• Simply out of curiosity it would be interesting to run the ML experiments again but on Colab Pro with a TPU connected, and see just how fast we could go.

Appendix

In this article we’ve explored how you can speed up training of ML models using the GPU in recent Mac devices via Apple’s Metal API. We’ve shown that this is possible with both TensorFlow and PyTorch. What we haven’t explored is how to actually enable GPU acceleration in either of these libraries! We’ll correct that here; detailing how to install each library with Metal integration enabled, and also how to utilise Metal acceleration in both TensorFlow and PyTorch code.

TensorFlow

There is a specific build of TensorFlow with Metal acceleration available from Apple, along with install instructions, which must be followed exactly. This is particularly true on the Intel Mac, where in testing I found that the build of TensorFlow with Metal acceleration was only found by pip when installed in a Python virtualenv — and not in a conda env. The build is only available for Python 3.8.

To check Metal acceleration is available to TensorFlow:

import tensorflow as tf
tf.config.list_physical_devices()

If you get the following result (or similar, but with a GPU in the list of PhysicalDevices), Metal support has been properly configured for TensorFlow:

[PhysicalDevice(name='/physical_device:CPU:0', device_type='CPU'),
 PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]

Note that you do not need to explicitly target the GPU for your workloads— TensorFlow will use it automatically if it is present.

PyTorch

It’s a lot easier to install PyTorch with Metal acceleration enabled — you just need to install v1.12 or later. At time of writing this requires installing PyTorch from the preview channel, but once PyTorch v1.12 moves out of pre-release this will get simpler again. Details are available in the PyTorch install instructions.

Conversely, it’s a little more difficult to check Metal acceleration is available in PyTorch applications and then also use it to accelerate training. You can use the following to check if Metal acceleration is available:

import torch

mps_built = torch.backends.mps.is_built()
mps_avail = torch.backends.mps.is_available()
device = "mps" if mps_built & mps_avail else "cpu"

Note that PyTorch refers to Metal acceleration as MPS (Metal Performance Shaders). To utilise the MPS device to train a model you can then do something similar to the following, having predefined model:

model = model.to(device)

Accelerated Data Science was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

We are happy to announce that our collaboration with Google DeepMind on applying deep learning techniques to precipitation nowcasting has now been published in Nature.

Skilful precipitation nowcasting using deep generative models of radar

This work has been made possible by leveraging a broad range of expertise from across the Met Office, from Operational Meteorology through Data Engineering, from User Experience researchers to Earth scientists, statistical scientists and many more. This diverse team was paired with Deep Learning and Machine Learning expertise from Google DeepMind to create a highly multidisciplinary team to deliver this groundbreaking research. The paper builds on our earlier review paper in which we summarise the emerging field combining Machine Learning and precipitation nowcasting.

A review of radar-based nowcasting of precipitation and applicable machine learning techniques (arXiv)

We are excited to have been able to contribute to this emerging field, and look forward to continuing to drive forward Environmental Intelligence research, exploring the intersection between machine learning and the Earth sciences.

Deep generative nowcasting: Nature publication was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

Climate scientists often use causal reasoning informally, but formalising it can have a couple of advantages. First, it makes it easier for others to grasp their assumptions and follow their argument. Second, it helps to apply the same reasoning to complex scenarios fluently and with ease, simplifying the analysis.

These advantages can be gained by moving from reasoning informally with causal narratives, to reasoning formally with causal diagrams.

Diagrammatic reasoning makes it possible to capture complex chains of reasoning by the application of a few simple visual rules. Once a climate scientist has developed a theory of the causal relationships of the indices they are considering, by drawing the network explicitly they can use these rules to determine which covariates to include and to exclude for their specific analysis. These rules are demonstrated by using examples from climate science, in the “Quantifying causal pathways of teleconnections” paper, which was recently published by BAMS (Bulletin of the American Meteorological Society).

The following researchers from Met Office’s Informatics Lab and University of Reading came together to collaborate on this paper:

Marlene Kretschmer (Department of Meteorology, University of Reading),

Samantha V. Adams (Met Office Informatics Lab),

Alberto Arribas (Microsoft Europe),

Rachel Prudden (Met Office Informatics Lab),

Niall Robinson (Met Office),

Elena Saggioro (Department of Mathematics and Statistics, University of Reading),

Theodore G. Shepherd (Department of Meteorology, University of Reading)

To highlight the major points of this paper, let’s take a look at the abstract:

“Teleconnections are sources of predictability for regional weather and climate but the relative contributions of different teleconnections to regional anomalies are usually not understood. While physical knowledge about the involved mechanisms is often available, how to quantify a particular causal pathway from data is usually unclear. Here we argue for adopting a causal inference-based framework in the statistical analysis of teleconnections to overcome this challenge. A causal approach requires explicitly including expert knowledge in the statistical analysis, which allows one to draw quantitative conclusions. We illustrate some of the key concepts of this theory with concrete examples of well-known atmospheric teleconnections. We further discuss the particular challenges and advantages these imply for climate science and argue that a systematic causal approach to statistical inference should become standard practice in the study of teleconnections.”

Purpose

In a nutshell, the purpose of this paper is not to detect causal relationships, but to make climate scientists aware of the need to bring explicit causal reasoning into their statistical analysis of teleconnections. Existing knowledge about physical mechanisms can be used to quantify teleconnection pathways and distinguish between direct and indirect pathways of influence. The paper includes some concrete examples of prominent teleconnections such as the influence of ENSO on California precipitation.

Summary

Throughout this paper, causal networks are used to represent teleconnections to facilitate quantification of the strength of the pathways. The authors advocate for a formal causal framework in the statistical analysis of teleconnections, which can be obtained by grounding it in causal inference theory (as seen in the paper’s info box). Statistical analysis in weather and climate science is usually done in the context of physical, hence causal, reasoning, but this reasoning is often only informal.

This is illustrated with a number of well-known teleconnection examples. Most of the examples are shown using Multiple Linear Regression (MLR), but all the concepts extend naturally into the non-linear context, as is demonstrated with a final example.

Some practical challenges for the use of a causal framework in climate science are also discussed in the paper by the authors. For instance, the various temporal scales of dependencies in the climate system can be difficult to address, as loops and cycles are generally not permitted in a causal network.

Another challenge is that relevant processes should be included in the network, meaning that no confounders exist there. This can never be fulfilled for an object as complex as the climate system. On the other hand, the same challenge (and limitation) exists for any statistical analysis.

Bearing in mind the concerns above, as regards to teleconnections, there is usually enough expertise available to draw a plausible causal model to articulate physical hypotheses. The authors argue that every climate science paper invokes, in some way, physical hypotheses, as any statistical study design depends on causal assumptions. Causal networks turn these claims into testable objects, making research conclusions more transparent and traceable.

Wrapping Up

The transparent and rational nature of causal network analysis can help in overcoming many of the limitations faced in current studies, and in reconciling differences between the conclusions of different studies. However, a causal approach is not meant to compete with traditional climate model experiments or physical theory. Instead it serves as a “scaffold” to build scientific intuition into the statistical analysis of the data. Both basic physics as well as data science are needed to make progress in climate science, and causal theory is a framework for better reconciling the two.

Many thanks to Rachel Prudden, Sam Adams, Marlene Kretschmer and Ted Shepherd for providing background information.

Further information and links:

• Quantifying causal pathways of teleconnections
• Bulletin of the American Meteorological Society
• Supporting material

Using causal network theory as a scaffold for scientific intuition was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

The 2021 Summer Meeting of the Earth Science Information Partners (ESIP) took place (virtually of course) from 19 to 23 July 2021. ESIP is a group of Earth Science organisations, currently primarily from the United States, that aims to promote “the collection, stewardship and use of Earth science data, information and knowledge that is responsive to societal needs” (from ESIP website). The group is organised around clusters, which are small interest or working groups, each with a particular focus within the broad remit of Earth Sciences. Some examples include the Cloud Computing cluster, which focus on the application on cloud-based- architectures to Earth Sciences, the Machine Learning cluster, focusing on best practices in applying Machine Learning to Earth Sciences, and the Data Readiness cluster, which is developing standards for evaluating how easily a dataset can be used for Machine Learning or Artificial Intelligence applications.

The clusters comes together twice a year for the Summer and Winter meetings to share knowledge and ideas between clusters as well as with the wider community. In this blog post I will take a look at some of the key themes and take away messages from my perspective from the sessions I attended and discussions I participated in, and consider what this means for the application Data Science and Machine Learning techniques in the Earth Sciences context.

At each conference there are a few key themes that emerge, popping up in almost every talk, representing the common challenges and opportunities in that field at that time. While distilling my conference experience into a few key themes, the words of the song “Nowadays” from Chicago popped into my head:

There’s men
Everywhere jazz, everywhere booze
Everywhere life, everywhere joy
Everywhere, nowadays

Lyrics of Nowadays from musical Chicago composed by John Kander & Fred Ebb (DOI not available)

The musical Chicago represents a slightly glamourised, slightly cynical view of life and fame in the roaring twenties, a golden age in many ways before the Great Depression. We’re currently in a golden age in the adoption of data-driven techniques. The talks and discussions at the ESIP 2021 Summer Meeting represent a snapshot of the state of Data Science and Data engineering in the Earth Sciences domain, both where we are now and where we would like to be in 2021. What key technologies and principles are “everywhere”, at least on the minds of everyone attending the ESIP meeting? What does the Data Science for Earth Sciences landscape look like “nowadays”?

Data everywhere

Sessions:

• OpenDAP for Data Providers (recording)
• AI Data Readiness: Designing A Community-Driven Road Map for Data Standards and Tools (recording)
• The Saga Continues: Cloud-Optimized Data Formats (recording)

Ideas & Tools:

• Streaming data, optimising for cloud access e.g. OpenDAP, EDR
• Making data easier to find through catalogs and search e.g STAC

It’s hardly novel for data to be a key theme of an Environmental Sciences conference, as data from both observations and simulations have long been essential. New data sources though do continue to be developed and existing data sources are made easier to access. In the new sources categories, we have sensors like “smart webcams”, which use Machine Learning to infer measurements like cloud cover from images or video. There has also been a lot of effort devoted to finding data sources through data catalogs, using technologies like STAC or EDR, as well as data search engines and portals, such as NASA GIBS or Radiant ML Hub. We see increasingly in applications that value from the large quantity and variety of data produced not from its use in isolation for a particular purpose, but used together with many different datasets from different domains, in ways the producer may not have expected. These users may often then not be experts in the domain of the data and data producer. Which brings us on to the next key theme: users.

Users everywhere

Sessions:

• Delivering Trusted Data to Real Users & Decision Makers (recording)
• Distributed Rapid Collaboration on Disaster related Information (recording)
• New Frontiers in AI for Earth and Space: Big Data and Parallel Computing (recording)

Ideas & Tools:

• Data Readiness levels & Decision-ready Information
• Different Data Use Cases e.g. GIS user, disaster responders

What do users really require from our data? Many users are joining together data from different sources and domains to draw insights, so it can be difficult for data providers to understand all user. This doesn’t mean as data producers or curators we should throw up our hands and give up trying to understand our users or how the data is used. Rather it is an ongoing key task to monitor who is using the data, for what and how, and letting that inform the contents of datasets and processing tools that use the data so that they meet actual users’ needs as much as possible.

For example, in the Analysis-ready Data cluster I am a part of, we are putting together a checklist about how ready a dataset is to be used for AI applications. Although some properties apply to all datasets and users, we need to consider specific use cases and the usefulness of our checklist for specific users to create a more useful document.

One user category that is fairly new to me as a Data Scientist is thinking about users like disaster responders, who don’t want access to huge amounts of raw data, but rather want to be presented with decision-ready information. They need to know though to what degree they can trust this information and whether they can have confidence in the decisions they make based on that information? How do we present our datasets and build pipelines so that users can trust the information products that come out of the end?

FAIR everywhere

Sessions:

• Best Practices for FAIR Research Software (recording)
• Best Practices for Reusability of Machine Learning Models: Guideline and Specification (recording)

Ideas & Tools:

• Apply FAIR to all assets
• Catalog tools — STAC, EDR
• Data discovery / search engines — e.g. NASA GIBS, Radiant ML Hub

Many readers will have realised that the challenges described so far are addressed by the FAIR principles — Findable Accessible Interoperable Reusable. I’m familiar with FAIR applied to data. What was new to me at this conference was extending the application of these principles to all digital assets in the processing pipeline, including the code and machine learning models.

Making code or ML models findable and accessible is fairly straightforward. Code has git repositories or zenodo archives. Standard pathways are less defined for ML models and not at the level of data search engines and catalogs described earlier, but are gradually evolving and developing with services like DL Hub.

Interoperability is perhaps better defined for code and models. Code should have a well defined interface (API) so that it can be joined with other analysis code. Models should have well defined inputs and output to join with other models or analysis code. Reusable is a more contested concept. What sort of reuse do we expect? Do we really mean reproduce or replicate, terms which are often used imprecisely or interchangeably, but mean different things. Technologies like conda environments and docker containers help to make code and models reusable and often reproducible, especially in association with FAIR data. Formats now exist for exchanging and reusing trained model parameters such as ONNX.

Ultimately what the FAIR principles mean for each asset and context comes down to what the expected users and the defined use cases to be supported. One can never perfectly apply FAIR for unknown uses, but hopefully in considering broad usage patterns one can get close for most users.

Metadata everywhere

Sessions:

• Machine-Readable Descriptors for Heterogeneous Tabular Data (recording)
• Toward Improving Representation of Data Quality Information (recording)

Ideas & Tools:

• Metadata vital for FAIR and ARCO, builds trust in data and derived products.
• Metadata for tabular — CSV for Web, CSV YAML, ERDDAP

We will get the most value out of our datasets when they can easily be used by users from different backgrounds and domains. It essential for these diverse users, who may not be from the domain of the data provider, to be able to get relevant information about using the data. All datasets therefore need to have comprehensive and comprehensible metadata. This includes describing the fields, their units, the origin of the data, the processing that has been applied and any other information necessary to understanding and using the data. This is also important for efficient access patterns, as tools which extract a subset or stream data need to know the detailed structure of the data from the metadata, without reading and of the actual data.

The challenge for providers is in finding the time for the labour intensive task of gathering and documenting the metadata along with all the other tasks required of data providers. Providers are also often primary data users, but often don’t need the metadata themselves as they know the data intimately, so curating metadata is devoting time to a task that will primarily benefit others, or possibly one’s future self. Metadata builds confidence for downstream products that they are built on correct use of datasets, and therefore confidence in the decisions made based on the derived decision-ready information.

Another challenge is that it is not always clear what is the best way to store the metadata so it is attached to the data. For example CSV files are a common way of storing tabular data, but do not have good innate support for detailed metadata. The session on Machine-Readable Descriptors for Heterogeneous Tabular Data suggested several emerging standards that extend the CSV format in a backward compatible way, such as CSV for the Web . It is important in building a solid foundation for a data driven future that metadata is not an afterthought or optional extra, but seen as vital and used extensively by data provider and user communities. Such tools help, but ultimately it is a culture shift that communities of interest and practice must foster, which leads onto the final theme of Community Everywhere.

Community Everywhere

Sessions:

• AI Data Readiness: Designing A Community-Driven Road Map for Data Standards and Tools (recording)
• ESIP in 2031: how we got here, from a pandemic to a bright new future (recording)

Key idea / tools / technologies:

• Community driven & owned recipes & standards
• Provide mentorship and training to increase uptake of these ideas and tools.

There were a lot of great ideas and tools in every session of the conference. I’ve never gone away from a conference with so many different things I wanted to try and so many tools that would be directly relevant to projects I’m currently working on. I also really enjoyed the space allowed for discussion, breakout sessions and socials throughout the conference. The organisers did an excellent job creating a community feel that encouraged real engagement.

After the post-conference high, though, you get back to the daily reality of one’s job with deadlines and deliverables and the high-minded ideals you intedned to adopt and the new tools you were eager to try are moved the bottom of your to-do list. To ensure that principles such as FAIR and tools such as metadata-rich cloud optimised formats and searchable data catalogs become a part of the usual working practicesin Earth Sciences, there needs to be a community to facilitate this and drive the process forward, developing community recipes and standards that make it easy to adopt best practices, along with mentorship and learning opportunities to develop the skills required.

One very interesting session was the last of the conference, which was a time-travel thought experiment. Imagine it is the 2031 and the ESIP community has grown and thrived over the past ten years of the “post-pandemic new normal”. What decisions were made and actions taken from the 2021 Summer Meeting and subsequently to make that happen? I think all conferences should include this combination retrospective / future planning exercise as a satisfying way to draw things to closure while ensuring the right “next steps” are taken. One overwhelming theme in this discussion was the importance of inclusive, active communities in making this happen.

Conclusions: The Future of Data Science for Earth Sciences

You can like the life you’re living, you can live the like you’re like.

Chicago (again)

What do we want the application of Data Sciences to Earth Sciences to look like in 10 years time? There is a lot of excitement about how technologies such as cloud computing and machine learning can revolutionize the domains of Earth Sciences. We’re moving beyond the initial excitement phase of “Look, I applied machine learning to my subject” or “I ran my script on the cloud”. We need to consider what can be improved or even transformed for real uses with specific use cases, and what principles and tools can provide a framework for turning those possible benefits into a reality.

The musical Chicago talks about the roaring twenties, an economic boom period, but the excesses of that era led to the devastating Great Depression that follows. Machine Learning, like most technologies, has gone through boom and bust cycles of hype and scepticism. We are currently going through a boom in adoption and excitement, and we want to ensure that through the maturation of these fields these themes represent that we don’t over promise and under deliver leading to a bust, but rather create real value to our organisations, communities and the wider world and deliver some of the transformational capabilities long promised and expected.

Further information and links:

• ESIP Home
• Clusters
• ESIP Machine Learning cluster
• ESIP Data Readiness cluster
• ESIP Data Lifecycle cluster
• ESIP Cloud Computing

ESIP Summer meeting 2021:

• Agenda
• Recordings of session — YouTube

The State of Data Science for Earth Sciences was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

Do give it a watch! I would also highly recommend browsing the Oxford ML and Physics Seminars YouTube channel, which is a treasure trove of fascinating talks on a huge variety of topics, from generative modelling of quantum states to seasonal sea ice forecasting.

The full abstract for my talk:

In both atmospheric science and machine learning, it is important to capture the uncertainty of predictions. This information can avoid the dangers of relying on over-confident predictions which may be incorrect, and help to understand the potential for high-impact rare events. Nonetheless, to focus only on capturing uncertainty risks giving an incomplete picture of the strengths of probabilistic modelling. At their heart, probabilistic models are about information. How does a new observation influence our mental model of the world? In this talk, I will discuss the use of probabilistic models in atmospheric science and how they can aid scientists in interpreting incomplete observations. I will also present some early work in this direction for spatial reasoning and time-series analysis.

Probabilistic Reasoning (Seminar) was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

As we’ll see, things get even more interesting when we combine both types of observation.

Why might we want to do this? Let’s look at a motivating situation. Suppose you have an autonomous vehicle recording a stream of observations relating to the rainfall at its current location. These tell you a lot about the situation along the vehicle’s past trajectory, but not much about what it will encounter further on its journey. On the other hand, radar observations will give you some information for the entire journey, but with a resolution that’s much lower. The question then is, can we gain any further insight by combining these two sources of information so that the vehicle can improve its awareness of the conditions it is likely to encounter?

In fact, something more interesting is true. At least in theory, there is information we can glean from the combination of these two types of observation that we couldn’t get from either individually. What’s more, it emerges from the mathematics describing the situation in a beautiful and natural way.

Illustrative example

To illustrate the idea, in this post I will use a toy model combining a GRF model with synthetically generated data and perfect observations. This will let us get into the interesting mechanics of the problem without getting hung up on thorny practicalities like noisy measurements and the non-Gaussianity of rainfall. Of course, these are important practicalities that would need to be considered in a real application.

First, we’ll generate some synthetic data. Here’s a sample from a Gaussian random field*, generated using the tools from the previous post:

Here’s the same data, but with reduced resolution:

We can use the techniques from the previous post to get a distribution for the high resolution data (assuming for a moment that we already have a good estimate of the prior lengthscale). If we sample from this distribution, we get something like this:

Now, what happens if we also have a point observation? To make the difference stand out, we’ll use an observation that’s quite a bit higher than its surroundings. This will let us see how the observations propagate through the mathematics, and affect the final result.

Below, I’ve shown a sample from a model incorporating a single point observation in addition to the low resolution averages, compared against a sample from our original model using low resolution averages only. The point observation has a value of 2, and its location is indicated by the red arrow. It’s clear that the point observation has pulled up the value of the field, not only at that point but in the surrounding area.

Analysing the distribution

To get a deeper understanding of what’s going on, we’ll need to go beyond looking at samples and visualise the distribution as a whole. We can do this by looking at the mean field and the variance at each point.

Shown below are the mean fields for the models with and without the point observation, and on the right the difference between the two.

Now we can see there’s a bit more going on than was immediately apparent. The point observation hasn’t just pulled up the surrounding values, it’s actually pushed down the values in the rest of the grid box.

The reason this happens is the interplay between the two kinds of observation. We’ve observed a high value at the marked point, meaning that the values in surrounding neighbourhood are also likely to be high. But we’ve also observed a low-ish average value for the grid box as a whole. To compensate for the high values in one corner, the values in the rest of the box must be lower. This effect even spills over into neighbouring gridcells; the point observation induces a higher value in one corner, and lower values elsewhere.

What’s nice about this is that we didn’t have to build this reasoning into our model. The logic arises from the standard linear algebra calculations we used to condition the model on our observations. If you remember back to the previous post, the pattern of negative values in the posterior covariance matrix for the spatial model was another manifestation of the same effect.

Now let’s have a look at the variance of the model distribution at each point. The first thing to jump out is the checked pattern present in both models. This isn’t an artefact, but reflects the information available from the spatial averages: these are more informative near the centres of the grid boxes, less so on the boundaries between them, and least of all at the corners where four of them meet. That said, they do somewhat obscure the influence of the point observation.

The difference field on the right tells a clearer story. The biggest reduction in variance is in the area surrounding the point observation, as expected. What’s more interesting is the reduction in variance at some distance from this point, most noticeably in the opposite corner of the grid box. Comparing this pattern to the mean field diff above, the reduction in variance lines up closely with changes, positive or negative, in the mean field.

Wrapping up

Taking a step back, what does all this tell us? By combining two types of observation, a point value and a spatial average, we have gained information; this is shown by a reduction in variance. The neat part is, some of this information gain happens at some distance from the point observation, such as in the opposite corner of the grid box. This is information beyond what we’d already get from just the spatial observations. It is an emergent effect, in that neither the spatial nor the point observation alone tells us anything about this corner. It’s only when we combine the two that we learn something new.

That’s all for this two-part series on spatial reasoning in random fields!

Spatial Reasoning in Random Field Models: Part II was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

The Environmental Futures & Big Data Impact Lab is a partnership of seven organisations, including the Met Office and the University of Exeter. The Impact Lab is part-funded by the European Regional Development Fund. It provides support to Small to Medium Enterprises (SMEs) in South West England to help them make better use of big data tools and environmental data to grow their business. The Informatics Lab as part of the Met Office contributes Data Science and Weather and Climate expertise to Impact Lab projects and events.

One of the ways the Impact Lab helps SMEs to grow is by putting on events on particular topics, such as machine learning or utilising climate datasets, which may be pertinent to their business, but which they lack necessary expertise to start exploiting. These events can lead to SMEs approaching the Impact Lab with a project proposal to help them in a particular area.

In pre-pandemic times, these events would be held at the Exeter Science Park, where the Met Office and University of Exeter staff members of the Impact Lab are based, and would present an opportunity for networking. Remote working has changed the nature of these events, but also provided some advantages, such as the ease with which these events can be recorded for later viewing by those who could not attend the live presentation. In this blog post I will describe two of the recent virtual events that we have held. There are links to a recording of the event on YouTube and example code on GitHub. Information on upcoming events can found on the Impact Lab Eventbrite page.

Introduction to Python for Data Analysis

Following on from previous events which gave high level overviews of Machine Learning and Computer Vision, feedback from participants suggested that they were interest in applying techniques from these fields, but did know where to get started. In particular many said they had some programming experience, but not in python which is the most commonly used programming language in these fields, which presented a barrier to entry. This suggested we should put on an event which showed participants how to get started using python for data analysis.

The event assumes some basic programming knowledge and preferably some basic level python knowledge. The tutorial recaps the basic concepts of python and how one might handle data using basic python language features. It then covers use of more sophisticated data science tools available in Python, such as Pandas and Iris, to write more efficient and readable code for data handling, looking at concepts like metadata and tabular vs gridded data.

Recordings:

Example Code: https://github.com/informatics-lab/intro_python_data_analysis

Using Climate Data for Business Planning

Climate Change will affect many aspects of life in the coming century. Whereas climate research was viewed as primarily of scientific interest when the Met Office Hadley Centre was founded 30 years ago, it is now recognised that all organisations, including businesses big and small, need to factor the impacts of climate change and climate variability into their long term business planning. To consider climate impacts, SMEs need to know how and where to access public climate data sets and what tools they should use to extract and process data relevant to their business.

This tutorial gives a brief introduction to climate change and the different types of climate data sets, including observations and simulation output. We consider how changes to meteorological variable such as temperature translate into real-world impacts and what sort of planning questions might factor in the impact of climate change . We then show some examples of how to do basic analysis of climate change impact using the UK Climate Projections data set for use cases such as energy production or the cooling and heating requirements of buildings might change.

Example Code: https://github.com/informatics-lab/UsingClimateData

Further Information

• Impact Lab website
• Upcoming Impact Lab Events (eventbrite)

Impact Lab Data Science Events was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

Most of the time, like any other ML method, the input data used to train¹ these models consists of a bunch of isolated examples. Having access to the function values at these input points, the challenge is then to predict the values at unobserved locations. So far, so familiar.

But these models are more flexible than they first appear. A part of my PhD research has been looking into what happens when these models are used with fundamentally different kinds of observations, such as integrals and spatial averages. This turns out to be a pleasingly natural extension, and what’s more it gives rise to some very beautiful (and useful!) emergent effects.

I’m going to try and illustrate a few of these effects in this short two-part series of blog posts. In this first post, I’ll give some background on the mathematics behind these models and show how they can be combined with spatial averages for certain kinds of downscaling / super-resolution. In Part II, we’ll see how combining point observations with spatial averages can make things even more interesting.

Gaussian processes

The intuition behind Gaussian processes is: they give you a distribution over the space of functions. Every sample you draw from the distribution is a function. The properties of the function, such as its smoothness, are encoded in the parameters of the distribution, through its covariance.

There are many excellent resources which explain the basics of Gaussian processes, so I won’t go into much detail here. If the ideas are new to you, I would definitely recommend watching an introductory talk like this one. Another great place to start is this article on distill, which walks through the mathematics with interactive examples.

As their name suggests, Gaussian processes are built on Gaussian distributions. In fact, if you take a finite² Gaussian process defined over a grid, it is just a multivariate Gaussian distribution, defined by a mean vector and a covariance matrix. What makes them special is the form of the covariance matrix, which reflects the spatial properties of the functions, including their dimension (see the image above). This is easiest to see for the one-dimensional case: the covariance peaks along the diagonal and then drops away, as each point is most closely correlated with its near neighbours and less so with those further away.³

The useful properties of Gaussian processes are also inherited from Gaussian distributions. The most important are:

Conditioning

Knowing the function values at certain points, we would like to get the conditional distribution over functions that match these observations. We can do this by conditioning the Gaussian process on our observations, which conveniently comes down to some matrix algebra:

for the conditional mean and covariance respectively, where t and o denote the target and observed locations, and obs is a vector containing the observed values.

Inference

Given a Gaussian process and an observed function, it’s useful to be able to calculate its likelihood as a function of the GP parameters. Taking gradients through this likelihood then allows us to infer the model parameters based on the observations. The (log) likelihood also reduces to matrix algebra:

where K is a parameterised covariance matrix, and x is the observed data.

Sampling

Finally, if we have a Gaussian process defined by its mean and covariance, we’d like to be able to treat it as a generative model and draw samples. This reduces to… you’ll never guess…

yes, more matrix algebra! Here w is a vector of white noise, and the half power denotes a matrix square root (or usually the Cholesky decomposition, for efficiency).

That’s about it for Gaussian processes. Before we move on, a bit of nomenclature. Gaussian processes which have an interpretation as spatial fields (instead of being defined over time or some more abstract space) are often called Gaussian random fields or GRFs. Mathematically, it’s really all the same.

Spatial conditioning

I mentioned above that GPs can be used with different kinds of observations, such as spatial averages. So, how does this work? Once again, it comes down to the nice properties of the Gaussian distribution, specifically the fact that the covariance is linear. This means it’s possible to derive the covariance of spatial averages from the covariance of the original field. What we end up with is a distribution relating the spatial averages to the original variables. Once we have this, it’s easy to condition on spatial averages using exactly the same techniques described in the previous section.

This is still sounding quite abstract, so let’s look at an example. We’ll stick to a single dimension for now, to make the interpretation easier. Our data just consists of a series of values, shown below. We’ll use this data to see the difference between conditioning on point values and on spatial averages by using it as input to two different models. The first model treats the data as point observations, while the second treats them as observations of spatial averages. In all other respects the models are identical; both use the same prior distribution with zero mean and unit variance, and the same lengthscale parameter.

When we condition the first model on the input data, treating the observations as point values, we get something like this:

The distribution is constrained to pass exactly through each of the observed points, so the variance drops to zero at these points. Midway between observations the variance is much higher, and the distribution is more influenced by its prior. A similar picture arises from plotting samples from the distribution: each passes exactly through the marked points, but elsewhere they can be quite divergent.

How about our other model? In this case, there is no sharp decrease in variance since the observations are no longer located at a point. Instead, the variance stays fairly close to uniform throughout.

Similarly, if we look at samples from this model it’s clear that they are no longer constrained to go through the points marked with red circles, as in the previous model. Some samples go through them by chance, while others miss entirely.

Since this model does not constrain the samples to pass through some specified points, does it constrain some other property? The answer is yes, but it’s less easy to see from these plots. To see it more easily, let’s take a look at the spatial averages of the samples from each model. The samples from the model using point observations are shown in pale green, and the samples from the model using spatial observations are shown in dark green:

Now it’s clear that the spatial averages of samples from the second model agree exactly with each other, and with the input observations. By contrast, the samples from the first model have quite divergent spatial averages.

It’s clear that how we choose to interpret this input data can have a significant impact on our resulting distribution. But why? We saw above that conditioning a Gaussian process comes down to just some linear algebra. How does this linear algebra force our models to satisfy these different constraints?

One way to gain some insight into this question is to examine the covariance matrices. An advantage of working in one dimension is that these matrices are quite easy to interpret: each axis of the covariance matrix corresponds to the single spatial dimension, with the variance for each point along the diagonal. By examining the posterior covariance for each model, we can see how they encode the constraints we have imposed.

Consider the arrangement of positive, negative, and zero covariance values in each case. One notable feature is the presence of zero values on the diagonal of the point-based model, corresponding to the observed points. By contrast, the spatial model has positive values everywhere on its diagonal. This matches up neatly with our previous discussion.

Other features show up in the covariance matrices that would be difficult or impossible to see in our earlier visualisations. One example is the appearance of a positive-negative “dipole” around each observation point in the covariance of the point-based model. This feature can be seen as imposing smoothness of the function through the fixed point: if the function is higher on one side, it must be lower on the other to maintain continuity.

For the spatial model, the negative values in the covariance follow a noticeably different pattern. In fact, they follow closely the boundaries of the blocks over which we have taken the spatial averages. In this case the meaning is not to impose continuity, but to constrain the spatial average over these areas to match the observation. If there is a high value in one part of a block, the function elsewhere in that block must take on lower values to compensate, and vice versa. This is a theme we’ll pick up again in Part II of this series.

Wrapping up

That’s it for this post. We’ve gone through some background on Gaussian processes, and seen how conditioning them on spatial averages instead of point observations can change the properties of our model. In the next post, I’ll discuss combining spatial averages with point observations and how this leads to some nice emergent properties.

If you’d like to learn more about the ideas in this post there is now a preprint available which goes into greater detail and includes an application to super-resolution of weather fields. You can also watch my recent talk.

Notes

[1] In fact, the normal distinction between training a model and zero-shot learning gets a bit tricky here. If the input data to a GP represents house prices, it looks a lot like training a regression model; if the data represents pixel values, it has more of a zero-shot learning flavour.

[2] A nice feature of Gaussian processes is that you can also do the same things with continuous functions and ungridded observations, but we don’t need that here.

[3] The covariance matrix is usually taken to follow a simple parameterised form known as a kernel — this simplifies model inference, as only a small number of parameters need to be estimated.

Spatial Reasoning in Random Field Models: Part I was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.

How next-gen data platforms stop weather and climate scientists from being software engineers and focus on science

I recently gave a talk on the importance of data engineering, how it can make scientists more productive and introduced some exciting projects in the data engineering space. Please have a watch below, it’s only 15 minutes.

The talk was originally at IS-ENES3/ESiWACE2 virtual workshop on New Opportunities for ML and AI in Weather and Climate Modelling (16–18 March 2021). The IS-ENES3 YouTube channel has recordings of all the talks including an excellent one from my colleague Rachel Prudden on Stochastic downscaling using Gaussian random fields.

Any thoughts always appreciated below.

You do you was originally published in Met Office Informatics Lab on Medium, where people are continuing the conversation by highlighting and responding to this story.