hydroTSM (version 0.3-5)

hydrokrige: 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, ...)
## S3 method for class 'default':
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, 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.,

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 locat

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

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 <

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 co

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

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 i

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 t

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:

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 maxdi

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

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

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', '

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. Valid values are: -) R : uses the default setting of spplot -) numeric: vector of reals giving the exact values in which

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 -) numeric: a single numeric value that will replace all the NA values in x

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, it indicates the index of the co

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

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 packa

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.

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.

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://r-spatial.sourceforge.net/gallery/

See Also

krige, autoKrige, readShapePoly, spsample

Examples

Run this code

############
## 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 
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")
           
##################################################
## 2) Ordinary Kriging interpolation and plot, in catchments defined by a shapefile
## Probably you will need to resize your window
# 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)

##################################################
## 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
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)

##################################################
## 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

# Computing KED
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)

##################################################
## 5) Block IDW interpolation and plot of 'EbroPPtsMonthly' for 3 months
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")

Run the code above in your browser using DataLab