Alexander Barth, Aida Alvera Azcárate, Pascal Joassin, Jean-Marie Beckers, Charles Troupin.
Available under the Creative Commons Attribution-ShareAlike 4.0 license. $\newcommand{\vect}[1] {{\mathbf #1}}$
Gridding is the determination of a field $\phi(\vect{r})$, on a regular grid of positions $\vect{r}$ based on arbitrarily located observations. Often the vector $\vect{r}$ is on a 2D, 3D or even 4D space.
$\rightarrow$
This notebook uses Julia 1.0 and versions above, which can be downloaded from https://julialang.org/downloads/, and the following packages:
To install these packages, enter the Pkg REPL-mode (hit ] in the julia prompt) and type
add NCDatasets PyPlot IJulia DIVAnd
using NCDatasets
using PyPlot
using Statistics
Stored in the files dan_field_obs.nc
are the following variables:
f
: the (full) fieldF
: the observations (without noise)Fe
: the observations (with noise)X
, Y
: the location of the observationsx
, y
: the grid point of the field f
mask
: land-sea mask (true corresponds to sea and false to land)pm
: the inverse of the resolution in the x
directionpn
: the inverse of the resolution in the y
directiondatafile1 = "../data/dan_field_obs.nc"
if isfile(datafile1)
@info("File already downloaded")
else
download("https://dox.uliege.be/index.php/s/96B8MOQeIcaRUoV/download", datafile1)
end
[ Info: File already downloaded
ds = Dataset(datafile1)
# load field
f = ds["f"].var[:,:]
x = ds["x"].var[:,:]
y = ds["y"].var[:,:]
mask = ds["mask"].var[:,:] .== 1
# load observations
F = ds["F"].var[:]
Fe = ds["Fe"].var[:]
X = ds["X"].var[:]
Y = ds["Y"].var[:]
close(ds);
pm = ones(size(mask))*5;
pn = ones(size(mask))*5;
@show size(x);
size(x) = (50, 50)
f2 = copy(f);
f2[.!mask] .= NaN;
pcolor(x,y,f2; cmap = "jet");
colorbar();
(Hypothetical) example of oceanographic field.
e.g.
a peninsula or a dike).This is a plotting function that we will reused later
# common colorbar range
caxis = [16.5,22]
function plotfield(f2; ca = caxis)
f2 = copy(f2);
f2[.!mask] .= NaN;
contourf(x,y,f2,20; cmap="jet", vmin = ca[1], vmax = ca[2]);
#clim(ca);
colorbar()
plot(X,Y,"k.");
return nothing
end
# plot the field and compare to the true solution
function plotfieldcmp(f2)
figure(figsize=(15,6))
subplot(1,2,2)
plotfield(f)
title("true field")
subplot(1,2,1)
plotfield(f2)
end
plotfieldcmp (generic function with 1 method)
plotfield(f);
title("Field with the location of the observations");
Values of the true field extracted at the location of the observations
contourf(x,y,mask; levels = [0,0.5], colors = [[.5,.5,.5]])
scatter(X,Y,20,F; cmap="jet")
clim(caxis)
colorbar();
structures and fronts is lost
more information about its structure we include in the analysis, the closer we can get to the true field.
Observations are in general affected by different error sources and other problems that need to be taken into account:
Quality control is an important step to exclude suspicious data from the analysis. But since this is a subjective decision, the data should never be deleted but flagged as suspicious or bad data.
Let's plot the observations with noise (Fe
):
contourf(x, y, mask; levels = [0,0.5], colors = [[.5,.5,.5]])
scatter(X, Y, 20, Fe; cmap = "jet")
clim(caxis)
colorbar();
Because observations have errors, it is always better to produce a field approximation and never a strict interpolation.
Gridded field using linear interpolation. This method is
implemented in the function griddata
of MATLAB and GNU Octave. Right panel shows true field.
This figure shows what would happen if the observations would have been interpolated linearly.
The domain is decomposed into triangles where the vertices are the location of the data points based on the Delaunay triangulation.
Within each triangle, the value is interpolated linearly.
The noise on the observations is in general amplified when higher order schemes, such as cubic interpolation, are used.
There are a number of methods that have been traditionally used to solve the gridding problem:
Isohyet (lines of constant precipitation) drawn by hand (from http://www.srh.noaa.gov/hgx/hurricanes/1970s.htm)
As opposed to the subjective method, objective analysis techniques aim to use mathematical formulations to infer the value of the field at unobserved locations based on the observations.
There are several ways to define the weights, which result in different gridding techniques. We will review in the following the most used gridding methods.
Cressman weights depend only on the distance between the location where the value of the field should be estimated and the location of the observation:
Cressman scheme: $$\begin{alignat}{2} w(d) &\;=\;& \frac{R^2 - d^2}{R^2 + d^2} \quad \mbox{if $d < R$} \\ &\;=\;& 0 \quad \mbox{if $d \ge R$} \end{alignat} $$
Barnes scheme: $$ w(d) = \exp\left( - \frac{d^2}{2 R^2} \right) $$
x_coord = -4:0.1:4
# distance to observation at x=0
d = abs.(x_coord)
R = 2;
w = (R.^2 .- d.^2) ./ (R.^2 .+ d.^2);
w[d .>= R] .= 0;
plot(x_coord,w,"b", label = "Cressman scheme")
R = 1;
w = exp.(-d.^2 ./ (2*R.^2));
plot(x_coord,w,"r"; label = "Barnes scheme")
legend();
The search radius is the typical control parameter and defines the length-scale over which an observation is used. This length scale can be made to vary in space depending on data coverage and/or physical scales. This parameter is chosen by the users based on their knowledge of the domain and the problem.
include("../cressman2.jl")
cressman2
R = 2.
fi = cressman2(R,X,Y,Fe,x,y);
plotfieldcmp(fi); title("Gridded field by Cressman weighting");
The Cressman weighting is a very simple and numerically quite efficient method. However, it suffers from some limitations which are apparent in the previous figure.
No estimate can be obtained at locations when no observation is located within the $R$.
In regions with very few observations, the method can return a discontinuous field.
The presence of barriers cannot be taken into account easily.
All observations are assumed to have a similar error variance since the weighting is based only on distance.
As a variant of the Cressman weights, other weighting functions can be defined. In the Barnes scheme, the weights are defined using a Gaussian function
fi = cressman2(R/2,X,Y,Fe,x,y; method = :barnes);
plotfieldcmp(fi); title("Gridded field using Barnes weights");
variable at different locations are related in average__.
Hypotheses concerning those errors:
The optimal interpolation (OI) scheme can be derived as the Best Linear Unbiased Estimator (BLUE) of the true state which has the following properties:
From these assumptions the optimal weights of the observations can be derived to grid the field. The optimal weights depend on the error covariances. In particular,
diva_intro_slides.pdf
if you want to know how the weights are computedinclude("../OptimInterp.jl");
R = 2.2 * 1.6855;
R = 3.;
Var = fill(1/3,(length(Fe),));
meanF = mean(Fe)
Fa = Fe .- meanF;
fi,vari = OptimInterp.optiminterp((X,Y),Fa,Var,(R,R),length(Fa),(x,y));
fi = fi .+ meanF;
As an illustration, the field obtained by optimal interpolation is shown
plotfieldcmp(fi); title("Optimal interpolation");
Relative error variance of the gridded field by optimal interpolation
plotfield(vari); clim(0,1)
The variational inverse method (VIM) or variational analysis (Brasseur et al., 1996) aims to find a field with the following characteristics:
This is achieved by defining a cost function which penalizes a field that doesn't satisfy those requirements.
Error fields can be derived since:
The tools DIVA and DIVAnd implement the variational inverse method.
Field reconstructed using DIVA implementing the variational inverse method. Right panel shows true field.
using DIVAnd
┌ Warning: The active manifest file has dependencies that were resolved with a different julia version (1.9.2). Unexpected behavior may occur. └ @ ~/Projects/Diva-Workshops/Manifest.toml:0 Precompiling DIVAnd ✓ BitFlags ✓ Adapt ✓ LogExpFunctions ✓ LoggingExtras ✓ URIs ✓ ConcurrentUtilities ✓ Tables ✓ WoodburyMatrices ✓ TranscodingStreams ✓ OpenSpecFun_jll ✓ AlgebraicMultigrid ✓ EzXML ✓ Mustache ✓ AxisAlgorithms ✓ CodecZlib ✓ RecipesBase ✓ ChainRulesCore ✓ OpenSSL ✓ IterativeSolvers ✓ LogExpFunctions → LogExpFunctionsChainRulesCoreExt ✓ HTTP ✓ StaticArrays ✓ SpecialFunctions ✓ Adapt → AdaptStaticArraysExt ✓ StaticArrays → StaticArraysStatisticsExt ✓ OffsetArrays ✓ SpecialFunctions → SpecialFunctionsChainRulesCoreExt ✓ Interpolations ✓ DIVAnd 29 dependencies successfully precompiled in 14 seconds. 54 already precompiled.
Type ?DIVAnd
for the documentation of DIVAndrun.
len = 2.
epsilon2 = 1/3;
fi,s = DIVAndrun(mask,(pm,pn),(x,y),(X,Y),Fa,(len,len),epsilon2);
fi = fi .+ meanF;
plotfieldcmp(fi); title("variational inverse method (DIVAnd)");
Relative error of the field reconstructed using DIVA
err = DIVAnd_cpme(mask,(pm,pn),(x,y),(X,Y),Fa,(len,len),epsilon2);
plotfield(err; ca = [0,1]);