hydrokrige {hydroTSM}R Documentation

Krige for Hydrological Time Series

Description

Automatic interpolation for hydrological ts, with optional plot (wrapper to some functions of the gstat and automap packages)

Originally it was thought as a way to make easier the computation of average precipitation over subcatchments (given as input in a shapefile map), based on values measured at several gauging stations, but nowadays it can also be used for interpolating any variable over a grid given by a raster map.

Available algorithms: inverse distance weighted (IDW), ordinary kriging (OK) and kriging with external drift (KED)

The (Block) Inverse Distance Weighted (IDW) interpolation is a wrapper to the idw function of the gstat package (so, it requires the gstat package).

The automatic kriging (OK or KED) is a wrapper to the autoKrige function of the automap package (so, it requires the automap and gstat packages), which automatically selects the best variogram model from four different ones: spherical, exponential, gaussian and Matern with M. Stein's parameterization (for more details, see autoKrige)

Usage

hydrokrige(x.ts, x.gis, ...)

## Default S3 method:
hydrokrige(x.ts, x.gis, X= "x", Y= "y", sname, bname, 
        elevation, predictors, catchment.name = "all", type="cells", 
        formula, subcatchments, IDvar = NULL, 
        p4s=CRS(as.character(NA)), cell.size = 1000, 
        grid.type = "regular", nmin = 0, nmax = Inf, maxdist = Inf, 
        ColorRamp = "PCPAnomaly", plot = TRUE, col.nintv = 10, 
        col.at = "auto", main, stations.plot = FALSE, stations.offset,
        arrow.plot = FALSE, arrow.offset, arrow.scale, 
        scalebar.plot = FALSE, sb.offset, sb.scale, verbose = TRUE,
        allNA.action="error", ...)

## S3 method for class 'data.frame'
hydrokrige(x.ts, x.gis, X= "x", Y= "y", sname, bname, 
        elevation, predictors, catchment.name = "all", type = "block", 
        formula, subcatchments, IDvar= NULL, 
        p4s=CRS(as.character(NA)), cell.size = 1000, 
        grid.type = "regular", nmin = 0, nmax = Inf, maxdist = Inf, 
        ColorRamp = "PCPAnomaly", plot = FALSE, col.nintv = 10,
        col.at = "auto", main, stations.plot = FALSE, stations.offset,
        arrow.plot = FALSE, arrow.offset, arrow.scale, 
        scalebar.plot = FALSE, sb.offset, sb.scale, 
        verbose = TRUE, allNA.action="error",
        dates=1, from, to, write2disk = TRUE, 
        out.fmt= "csv2",
        fname = paste(ColorRamp, "by_Subcatch.csv", sep = ""), ...)

Arguments

x.ts

numeric or data.frame object, with measured values at several locations.

-) x.ts CAN contain as many points as you want, e.g., all the gauging stations in the your database.
2) x.ts HAVE to contain -at least- some points (e.g., gauging stations) that are also present in x.gis

The names of each value in x.ts are used as the ID of each measurement point. When x.ts is a vector, this can be checked with names(x), whereas when x.ts is a data.frame, it can be checked with colnames(x.ts). The IDs of each measurement point have to be equal to the ID stored in the field sname of x.gis.

When x.ts is a data.frame, its structure have to be as follow:
-) 1st column : OPTIONAL. It MAY contain the dates, date-time or IDs that identify the measured values of each row. See dates argument.
-) 2nd...Nth column: Measured values at each point (e.g., gauging station). The name of the columns is used as the ID of each station -starting with a letter!!-.

x.gis

data.frame with the spatial information for each measurement point (e.g., gauging stations).

-) x.gis MAY contain as many points as you want, e.g., all the stations in your database
-) x.gis HAVE to contain -at least- the location of those stations that will be used for the interpolations

The MINIMUM fields that have to be present in this file, and their corresponding column index are those described by: X, Y, sname.

The ID of each measurement point, given in the field sname, has to be equal to the corresponding ID used in x.ts

X

character, field name in x.gis that stores the easting coordinate of each measurement point.

Y

character, field name in x.gis that stores the northing coordinate of each measurement point.

sname

character, field name in x.gis that stores the ID of each measurement point. the name of each measurement point HAS to start by a letter!!-

bname

OPTIONAL. character, field name in x.gis that stores the name of the subcatchment in x.gis that will be analysed.
ONLY necessary when catchment.name is different from all.

elevation

OPTIONAL. character, field name in x.gis that stores the elevation of the gauging stations (m a.s.l.).

predictors

OPTIONAL. SpatialGridDataFrame-class object, with prediction/simulation locations (it is needed for KED).
Usually, a digital elevation model (DEM) read with the readGDAL function of the rgdal package. See the newdata argument in krige.
It should contain attribute columns with the independent variables (if present) and the coordinates with names as defined in x.gis
If predictors is missing, the grid to be used as prediction/simulation locations is generated from sampling the polygon specified by the user in subcatchments, according to the arguments provided by cells.size and grid.type

catchment.name

name of the catchment that will be analysed. Possible values are:
-)all : ALL the stations in the x.gis will be used
-)other string: ONLY those stations in x.gis with a bname field value equal to catchment.name will be used .

type

Character, indicating the type of interpolation required by the user.
When x.ts is a data.frame, the ONLY possible value is block. For all the other cases, possible values are:
-) cells : the interpolated values are computed for each cell
-) block : the interpolated values are computed for each subcatchment, where the value for each subcatchment is computed as the mean value over all the cells that belong to each subcatchment
-) both : cells and block are computed.

formula

OPTIONAL. Formula to be used in case of ordinary kriging or kriging with external drift. Requires the automap package. All the variables to be used within formula has to be present both in x.gis and predictors. See the formula argument in krige.
formula defines the dependent variable as a linear model of independent variables. Within this function, the dependent variable is always called value, therefore, for ordinary and simple kriging use the formula value~1; for simple kriging also define beta; for universal kriging, suppose value is linearly dependent on x and y, use the formula value~x+y.

subcatchments

OPTIONAL. Only required when predictors is missing.
Spatial polygon with all the subcatchments to be used as interpolation domain. The polygons are used to create the grid that will be used as prediction/simulation locations, from sampling it according to the arguments provided by cells.size and grid.type. Valid values are:
1) Character, indicating the filename (with path) of the shapefile that contains all the subcatchments to be used as interpolation domain. It HAS TO BE of 'polygon' type
2) SpatialPolygonsDataFrame-class resulting from reading the shapefile (e.g., with the command readShapePoly of the maptools package) that contains all the subcatchments to be used as interpolation domain.

IDvar

See readShapePoly. a character string with the name of a field in the subcatchments shapefile containing the ID values used to identify each one of the subcatchments. When type is block, the values stored in this field will be used for labelling the computed values in each one of the subcatchments, therefore, if you don't provide this value, it could be difficult to identify which computed value corresponds to which subcatchment, because the ID is assigned by the readShapePoly function.

p4s

OPTIONAL. a character with information about the projection of the GIS files, usually created by the CRS function of the sp package.

cell.size

OPTIONAL. Only required when predictors is missing. Size of the cells [m] to be used for sampling the polygons specified by the user in subcatchments, to create a grid to be used as prediction/simulation locations .

grid.type

OPTIONAL. Only required when predictors is missing. See spsample. Character, indicating the type of grid to be computed over the area defined by subcatchments. Valid values are:
-) regular : for regular (systematically aligned) sampling; Default option
-) random : for completely spatial random;
-) stratified : for stratified random (one single random location in each "cell"
-) nonaligned : for non-aligned systematic sampling (nx random y coordinates, ny random x coordinates);
-) hexagonal : for sampling on a hexagonal lattice;
-) clustered : for clustered sampling

nmin

OPTIONAL. See krige. For local interpolation: if the number of nearest observations within distance maxdist is less than nmin, a missing value will be generated; see maxdist. By default nmin=0.

nmax

OPTIONAL. See krige. For local interpolation: the number of nearest observations that should be used for a kriging prediction, where nearest is defined in terms of the space of the spatial locations. By default, all observations are used.

maxdist

OPTIONAL. See krige. For local interpolation: only observations within a distance of maxdist from the prediction location are used for prediction or simulation; if combined with nmax, both criteria apply. By default, all observations are used.

ColorRamp

Function defining the colour ramp to be used for plotting the maps OR character representing the colours to be used in the plot. In the latter case, valid values are in: c('Precipitation', 'Temperature', 'PCPAnomaly', 'PCPAnomaly2', 'TEMPAnomaly', 'TEMPAnomaly2', 'TEMPAnomaly3')

plot

Logical, indicating if the interpolated values have to be plotted or not

col.nintv

integer, number of colours that have to be used for plotting the interpolated values

col.at

Specify at which interpolated values colours change.[MJ m-2 day-1] or [mm d-1] at <- seq(min, max,length.out=col.nintv)
min <- floor(min(idw["var1.pred"]@data, na.rm=TRUE))
max <- ceiling(max(idw["var1.pred"]@data, na.rm=TRUE))

main

Character with the title to be used for the plot.

stations.plot

Logical, indicating if the gauging stations, defined by x.gis have to be plotted

stations.offset

2D list with the numeric coordinates in which the label with the amount of gauging stations have to be plotted. e.g., stations.offset = c(450000, 4600000).

arrow.plot

Logical, indicating if a North Arrow have to be plotted

arrow.offset

2D list with the numeric coordinates in which the north arrow has to be plotted. e.g., arrow.offset = c(690000,4760000)

arrow.scale

Scale (in the map units) to be used for plotting the north arrow, e.g., scale = 20000

scalebar.plot

Logical, indicating if a Scale Bar has to be plotted

sb.offset

2D list with the numeric coordinates in which the North Arrow has to be plotted. e.g., sb.offset = c(400000,4490000)

sb.scale

Scale (in the map units) to be used for plotting the Scale Bar, e.g., scale = 100000, means that the scale bar will have a length of 100km

verbose

logical; if TRUE, progress messages are printed

allNA.action

Action to be executed when all the values in x.ts are NA. Valid values are:
-) "error": it will produce an error message. Default option
-) a single numeric value that will replace all the NA values in x.ts, giving place to a map with a constant value. At your own risk !

dates

numeric, factor, or character object indicating how to obtain the ID (usually dates) that will be used to identify the interpolation carried out for each row of x.ts.
If dates is a single number (default), it indicates the index of the column in x.ts that stores the dates
If dates is a factor or character vector, its values will be used as ID for the interpolations carried out in each row of x.ts.

from

Character indicating the starting date for the values stored in all the files that will be read.

to

Character indicating the ending date for the values stored in all the files that will be read.

write2disk

Logical. Indicates if we want to write the output into a CSV file. Default value is TRUE

out.fmt

OPTIONAL, only needed when write2disk==TRUE.
Character indicating the type of csv file to be written with the results. Valid values are csv and csv2. For more information, see write.table

fname

OPTIONAL. Character with the filename of the output file. Only needed when write2disk=TRUE

...

further arguments passed to or from other methods. In particular, these further arguments are passed to the function idw (gstat package) OR autoKrige (automap package), depending on the value passed to formula (see 'details' below):
-) for IDW, the arguments are passed to : idw(formula, locations, newdata, nmin, nmax, maxdist, ...)
-) for OK, KED, the arguments are passed to: autoKrige(formula, input_data, new_data, nmin, nmax, maxdist, verbose=verbose, ...)

Details

The type of interpolation (IDW, OK, KED) is obtained from the argument formula:

-) When formula is missing, an IDW interpolation, by calling the idw function in the gstat package, with formula = value~1.
-) When formula = value~1, an OK interpolation, by calling the autoKrige function, with formula = value~1.
-) When formula = value~pred1 + pred2 + ..., a KED interpolation, by calling the autoKrige function, with the formula specified by the user.

When type is block or both, a block interpolation is carried out for each subcatchment defined by subcatchments, so far, computing the average value over all the cells belonging to each subcatchment.

The automatic kriging is carried out by using a variogram generated automatically with the autofitVariogram function of the automap package.

Value

Cells

When type is cells, the output object is a SpatialPixelsDataFrame-class, which slot 'data' has two variables: 'var1.pred' and 'var1.var' with the predictions and its variances, respectively


Block

When type is block, the output object is a SpatialPolygonsDataFrame-class, which slot 'data' has four variables: 'x', 'y' with the easting and northing coordinate of the centroid of the catchments specified by subcatchments , and 'var1.pred' and 'var1.var' with the predictions and its variances, respectively


list(Cells, Block)

When type is both, the resulting object is a list, with the two elements previously described.

Note

IMPORTANT: It is you responsibility to check the validity of the fitted variogram !!.

Author(s)

Mauricio Zambrano-Bigiarini, mzb.devel@gmail

References

N.A.C. Cressie, 1993, Statistics for Spatial Data, Wiley.

Applied Spatial Data Analysis with R. Series: Use R. Bivand, Roger S., Pebesma, Edzer J., Gomez-Rubio, Virgilio. 2008. ISBN: 978-0-387-78170-9

Pebesma, E.J., 2004. Multivariable geostatistics in S: the gstat package. Computers & Geosciences, 30: 683-691

http://www.gstat.org/

http://rspatial.r-forge.r-project.org/

See Also

krige, autoKrige, readShapePoly, spsample

Examples

############
## Loading the monthly time series of precipitation within the Ebro River basin.
data(EbroPPtsMonthly)

## Loading the gis data
data(EbroPPgis)

## Loading the shapefile with the subcatchments
data(EbroCatchmentsCHE)

## Projection for the Subcatchments file
# European Datum 50, Zone 30N
require(sp)
p4s <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")

## Selecting the first day of 'EbroPPtsMonthly' for all the stations.
# The first column of 'EbroPPtsMonthly' has the dates
x.ts <- as.numeric(EbroPPtsMonthly[1, 2:ncol(EbroPPtsMonthly)])

## Setting the name of the gauging stations
names(x.ts) <- colnames(EbroPPtsMonthly[1,2:ncol(EbroPPtsMonthly)])

##################################################
## 1) IDW interpolation and plot
## Probably you will need to resize your window 
## Not run: 
x.idw <- hydrokrige(x.ts= x.ts, x.gis=EbroPPgis, 
                    X="EAST_ED50", Y="NORTH_ED50", sname="ID", bname="CHE_BASIN_NAME",
                    type= "both",
                    subcatchments= EbroCatchmentsCHE,
                    cell.size= 3000, 
                    ColorRamp= "Precipitation",	
                    main= "IDW Precipitation on the Ebro")

## End(Not run)
           
##################################################
## 2) Ordinary Kriging interpolation and plot, in catchments defined by a shapefile
## Probably you will need to resize your window
## Not run:  

# Computing OK, over of 3000x3000m, sampled withinthe subcatchments defined by 'subcatchments'
x.ok <- hydrokrige(x.ts= x.ts, x.gis=EbroPPgis, 
                   X="EAST_ED50", Y="NORTH_ED50", sname="ID", bname="CHE_BASIN_NAME",
                   type= "both", formula=value~1,
                   subcatchments= EbroCatchmentsCHE,
                   p4s= p4s, 
                   cell.size= 3000, 
                   ColorRamp= "Precipitation", 
                   main= "OK Precipitation on the Ebro", arrow.plot= TRUE, 
                   arrow.offset= c(900000,4750000), arrow.scale= 20000,
                   scalebar.plot= TRUE, 
                   sb.offset= c(400000,4480000), sb.scale= 100000)


# Getting the interpolated values in each polygon 
# ('var1.pred' field in the output data.frame object)
x.ok.block <- slot(x.ok[["Block"]], "data")

# Getting the interpolated values in each cell (SpatialPixelsDataFrame object)
# (the pedicted values in each cell are stored in the 'var1.pred' field of the 
#  'data' slot)
x.ok.cells <- x.ok[["Cells"]]

# Plotting the interpolated values in each cell
spplot(x.ok.cells, "var1.pred")

## End(Not run)

##################################################
## 3)  Ordinary Kriging interpolation and plot, in an area defined by a raster map.
## The raster map may be any \link[sp]{SpatialGridDataFrame-class} object, read with 
## the \code{\link[rgdal]{readGDAL}} function of the \pkg{rgdal} package or similar. 
## Probably you will need to resize your window

#Loading the DEM
data(EbroDEM1000m)

#Giving a meaningful name to the predictor
EbroDEM1000m$ELEVATION <- EbroDEM1000m$band1

# Saving memory
EbroDEM1000m$band1 <- NULL

# Computing OK, over the spatial grid defined by the DEM
## Not run: 
x.ok <- hydrokrige(x.ts= x.ts, x.gis=EbroPPgis, 
                   X="EAST_ED50", Y="NORTH_ED50", sname="ID", 
                   formula=value~1,
                   p4s= p4s, 
                   predictors=EbroDEM1000m,
                   ColorRamp= "Precipitation", 
                   main= "OK Precipitation on the Ebro",
                   arrow.plot= TRUE, 
                   arrow.offset= c(900000,4750000), arrow.scale= 20000,
                   scalebar.plot= TRUE, 
                   sb.offset= c(400000,4480000), sb.scale= 100000)

## End(Not run)

##################################################
## 4) Kriging with External Drift interpolation and plot
## Probably you will need to resize your window

#Loading the DEM
data(EbroDEM1000m)

#Giving a meaningful name to the predictor
EbroDEM1000m$ELEVATION <- EbroDEM1000m$band1

# Saving memory
EbroDEM1000m$band1 <- NULL

# Forcing the projection of the DEM to be the same of the 'subcatchments' argument
# (just because I know in advance they are the same)
proj4string(EbroDEM1000m) <- p4s

# Computing KED
## Not run: 
x.ked <- hydrokrige(x.ts= x.ts, x.gis=EbroPPgis, 
                    X="EAST_ED50", Y="NORTH_ED50", sname="ID", 
                    bname="CHE_BASIN_NAME", elevation="ELEVATION",
                    type= "cells", 
                    formula=value~ELEVATION,
                    subcatchments= EbroCatchmentsCHE,
                    predictors=EbroDEM1000m,
                    cell.size= 3000, 
                    ColorRamp= "Precipitation", 
                    main= "KED Precipitation on the Ebro",
                    arrow.plot= TRUE, 
                    arrow.offset= c(900000,4750000), arrow.scale= 20000,
                    scalebar.plot= TRUE, 
                    sb.offset= c(400000,4480000), sb.scale= 100000)

## End(Not run)

##################################################
## 5) Block IDW interpolation and plot of 'EbroPPtsMonthly' for 3 months
## Not run: 
x.idw <- hydrokrige(x.ts= EbroPPtsMonthly, x.gis=EbroPPgis, 
                  X="EAST_ED50", Y="NORTH_ED50", sname="ID", 
                  bname="CHE_BASIN_NAME",
                  type= "cells", #'both'
                  subcatchments= EbroCatchmentsCHE,
                  cell.size= 3000,
                  ColorRamp= "Precipitation",
                  arrow.plot= TRUE, 
                  arrow.offset= c(900000,4750000), arrow.scale= 20000,
                  scalebar.plot= TRUE, 
                  sb.offset= c(400000,4480000), sb.scale= 100000,
                  dates=1, 
                  from="1942-01-01", to="1942-03-01")

## End(Not run)


[Package hydroTSM version 0.5-1 Index]