CAIRS is a framework to reconstruct rain fields by assimilating signals of fundamentally different rain sensors.
In particular, the integration characteristics of sensors are explicitly considered. For example, non-recording standard rain gauges integrate over time and deliver information such as the daily rainfall sums. The rain-induced attenuation of micro wave links (MWL) can be used to measure the path-integrated intensities - an example of a sensor with spatial integration.
Sensor signals with different scales (e.g. continuous, binary) can be assimilated. Furthermore, CAIRS is formulated continuously in time and space. This is helpful because it enables a natural consideration of signals with irregular time-intervals.
The mathematical model is described in Scheidegger and Rieckermann (2014). The basic functionality and application is explained in this tutorial.
Note, CAIRS is still under development and the interface may change.
CAIRS is a Julia package. The first step is to download and install Julia version 0.6 or newer.
CAIRS can then be installed with the Julia command Pkg.clone()
:
Pkg.clone("git://github.com/scheidan/CAIRS.jl.git")
After that, CAIRS behaves like a normal package. For example, it can
be updated with Pkg.update()
.
First, the package CAIRS must be loaded. For convinience, it is also
recommended to load the packages Dates
and Distributions
:
using CAIRS
using Distributions
using Base.Dates
Every sensor must be characterized. In the simplest case a sensor measures the rain intensity at a point. In this case (the logarithm of) the signal distribution must be defined conditioned on the intesity at this coordinate:
function log_p_gauge(S::Float64, R::Vector) # non-linear continuous rain gauge
mu = 0.1+R[1]^2.0 # Note, the signal and can be non-linearly
# related to the rain intensity.
sigma = 0.005
## log of normal density, p(S|R)
logpdf(Normal(mu, sigma), S) # doesn't have to be normal
end
sensor_gauge = Sensor(log_p_gauge)
For integrating sensors, also the integration domain must be
specified. For example, a micro wave link (MWL) with length 6
may be
defined as:
function log_p_MWL(S::Float64, I::Float64)
R_mean = I/6.0
sigma = 0.1
## log of normal density, p(S|I)
logpdf(Normal(R_mean, sigma), S)
end
sensor_MWL = Sensor(log_p_MWL, Coor(6, 0, 0)) # integrates along a path of length 6
The prior of the rain field is modeled as Gaussian process (GP). A GP is described by a mean and a covariance function.
This functions can be specified by the user. The mean function returns
the prior mean of the rain intensity at a given coordinate. It must
take a single argument of type Coor
. The covariance function must
return the covariance of the rain intensities at two given point, given
by two arguments of type Coor
. Note, it is not checked if the
provided function is a valid covariance function!
However, helpers to construct valid functions are provided. The functions
mean_constant()
and cov_exponential()
create a simple constant
mean, and a separable gamma-exponential covariance function. Only the
parameters must be provided:
mean_GP = mean_constant(mean=2.0)
cov_GP = cov_exponential(sigma=10.0, # standard deviation of GP
l_spatial=1.5, # spatial correlation length
l_temporal=Minute(1), # temporal correlation length
gamma=1.0) # exponent for smoothness in [0, 2]
Other types of covariance functions will be added in future.
The next step is to import the signals. Every signal must have an
attached sensor. Signals can be constructed with the function
Signal
or more conveniently with add_signal()
.
Currently add_signal()
expected that the signals of every sensor are
stored in a separate file. The file must contain two columns:
- Column 1: date and time
- Column 2: signal values
## path to example data that come with the CAIRS package
path1 = joinpath(Pkg.dir("CAIRS"), "example", "data", "Sensor1.csv")
path2 = joinpath(Pkg.dir("CAIRS"), "example", "data", "Sensor2.csv")
sig = Signal[] # create an empty array for Signals
add_signal!(sig, # add signal to vector 'sig'
path1, # file name
sensor_gauge, # sensor
Coor(5, 6), # coordinate of the sensor
date_format="d.m.yyyy HH:MM:SS",
delim=',') # delimitation character
add_signal!(sig, path2,
sensor_MWL, # MWL link
Coor(4.2, 2), # coordinate of one end point of the sensor
0.9, # rotation around the point defined above in [rad]
date_format="d.m.yyyy HH:MM:SS",
delim=',')
Information about a signal can be printed with show
, e.g. show(sig[1])
.
Writing the sensor positions in a file is useful for plotting:
sensor2csv(sig, "sensor_coor.csv")
The location for which a prediction of the rain intesity is desired must be
defined as an Array
or Vector
of coordinates. Coordinates are
defined with Coor(x, y, time)
. Time can be a number or a DateTime
object.
### create a simple grid
nn = 20
loc_pred = [Coor(i, j, time)
for i=linspace(0, 10, nn), j=linspace(0, 10, nn),
time=DateTime(2013, 11, 22, 13, 15, 00) : Minute(1): DateTime(2013, 11, 22, 13, 20, 00) ]
This produced a regular grid, but the point could also be irregularly distributed. Also, not only predictions for coordinates but also for intesities integrated over a domain can be made. Domains are defined by the function Domain
.
The assimilation of the signals and the computation of the predictions are done with predict
.
R_pred = predict(loc_pred, # vector or array with locations for predictions
sig, # vector of signals
mean_GP, # mean function of prior
cov_GP, # covariance function of prior
n_sample_calib = 20000, # number of iterations of the Gibbs sampler
burn_in = 5000, # number of removed samples (and length of adaptation)
n_sample_pred = 6000, # number of samples for predictions
delta = Second(90)) # consider all signals within time 'delta'
# from prediction points
Write a summary of the samples in a file that is used for visualization:
summary2csv(R_pred, "rain_field.csv")
One possibility to visualize the result is to use R. A simple
R-script to produce rain maps comes with CAIRE. It requires that R and
the R-libraries lattice
, latticeExtra
and tripack
are installed.
pathRscript = joinpath(Pkg.dir("CAIRS"), "R", "compute_rain_map.r")
run(`Rscript $pathRscript rain_field.csv sensor_coor.csv out.pdf`)
Note, here it is assumed that Rscript
is in PATH.
Scheidegger, A. and Rieckermann, J. (2014) "Bayesian assimilation of rainfall sensors with fundamentally different integration characteristics" in WRaH Proceedings, Washington, DC.