hydrokrige: Krige for Hydrological Time Series
In hydroTSM: Time Series Management, Analysis and Interpolation for Hydrological Modelling
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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://rspatial.r-forge.r-project.org/
See Also
krige, autoKrige, readShapePoly, spsample
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
############
## 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)
hydroTSM documentation built on March 13, 2020, 2:23 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note Author(s) References See Also Examples
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | 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. -) The names of each value in When |
x.gis |
data.frame with the spatial information for each measurement point (e.g., gauging stations). -) The MINIMUM fields that have to be present in this file, and their corresponding column index are those described by: The ID of each measurement point, given in the field |
X |
character, field name in |
Y |
character, field name in |
sname |
character, field name in |
bname |
OPTIONAL. character, field name in |
elevation |
OPTIONAL. character, field name in |
predictors |
OPTIONAL. SpatialGridDataFrame-class object, with prediction/simulation locations (it is needed for KED). |
catchment.name |
name of the catchment that will be analysed. Possible values are: |
type |
Character, indicating the type of interpolation required by the user. |
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 |
subcatchments |
OPTIONAL. Only required when |
IDvar |
See |
p4s |
OPTIONAL. a character with information about the projection of the GIS files, usually created by the |
cell.size |
OPTIONAL. Only required when |
grid.type |
OPTIONAL. Only required when |
nmin |
OPTIONAL. See |
nmax |
OPTIONAL. See |
maxdist |
OPTIONAL. See |
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: |
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]
|
main |
Character with the title to be used for the plot. |
stations.plot |
Logical, indicating if the gauging stations, defined by |
stations.offset |
2D list with the numeric coordinates in which the label with the amount of gauging stations have to be plotted. e.g., |
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.scale |
Scale (in the map units) to be used for plotting the north arrow, e.g., |
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.scale |
Scale (in the map units) to be used for plotting the Scale Bar, e.g., |
verbose |
logical; if TRUE, progress messages are printed |
allNA.action |
Action to be executed when all the values in |
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 |
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 |
fname |
OPTIONAL. Character with the filename of the output file. Only needed when |
... |
further arguments passed to or from other methods. In particular, these further arguments are passed to the function idw (gstat package) OR |
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 |
Block |
When |
list(Cells, Block) |
When |
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://rspatial.r-forge.r-project.org/
See Also
krige, autoKrige, readShapePoly, spsample
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 | ############
## 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)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.