Package 'fields'

Description

Adds an image to an existing plot. Simple arguments control the location and size.

Usage

add.image(xpos, ypos, z, adj.x = 0.5, adj.y = 0.5, 
image.width = 0.15, image.height = NULL, col = tim.colors(256), ...)

Arguments

See Also

image.plot, colorbar.plot, image, tim.colors

Examples

plot( 1:10, 1:10, type="n")
data( lennon)

add.image( 5,4,lennon, col=grey( (0:256)/256))
# reference lines 
xline( 5, col=2)
yline( 4,col=2) 

#
# add lennon right in the corner beyond the plotting region
# 

par(new=TRUE, plt=c(0,1,0,1), mar=c(0,0,0,0), usr=c(0,1,0,1))
add.image( 0,0, lennon, adj.x=0, adj.y=0)

Description

Adds arrows at specified points where the arrow lengths are scaled to fit on the plot in a reasonable manner. A classic use of this function is to depict a vector field. At each point (x,y) we have a vector with components (u,v). Like the arrows function this adds arrows to an existing plot.

Usage

arrow.plot(a1, a2, u = NA, v = NA, arrow.ex = 0.05,  
     xpd = TRUE, true.angle = FALSE, arrowfun=arrows,...)

Arguments

Details

This function is useful because (u,v) may be in very different scales from the locations (x,y). So some careful scaling is needed to plot the arrows. The only tricky thing about this function is whether you want the true angles on the plot. For overlaying a vector field on top of contours that are the streamlines true.angle should be false. In this case you want u and v to be scaled in the same way as the x and y variables. If the scaling is not the same then the arrows will not look like tangent vectors to the streamlines. An application where the absolute angles are meaningful might be the hands of a clock showing different times zones on a world map. Here true.angle=T is appropriate, the clock hands should preserve the right angles.

See Also

arrows

Examples

#
# 20 random directions at 20 random points

x<- runif( 20)
y<- runif( 20)
u<- rnorm( 20)
v<- rnorm( 20)
plot( x,y)
arrow.plot( x,y,u,v) # a default that is unattractive 

plot( x,y, type="n")
arrow.plot( x,y,u,v, arrow.ex=.2, length=.1, col='green', lwd=2) 
# thicker lines in green, smaller heads and longer tails. Note length, col and lwd are
# options that the arrows function itself knows about.

Description

Discretizes a set of 2-d locations to a grid and produces a image object with the z values in the right cells. For cells with more than one Z value the average is used.

Usage

as.image(Z, ind=NULL, grid=NULL, x=NULL,weights=rep(1, length(Z)),
 na.rm=FALSE, nx=64, ny=64, boundary.grid=FALSE,  nrow=NULL, ncol=NULL,
 FUN = NULL)

Arguments

Details

The discretization is straightforward once the grid is determined. If two or more Z values have locations in the same cell the weighted average value is taken as the value. The weights component that is returned can be used to account for means that have different numbers (or precisions) of observations contributing to the grid point averages. The default weights are taken to be one for each observation. See the source code to modify this to get more information about coincident locations. (See the call to fast.1way)

Value

An list in image format with a few more components. Components x and y are the grid values , z is a nrow X ncol matrix with the Z values. NA's are placed at cell locations where Z data has not been supplied. Component ind is a 2 column matrix with subscripts for the locations of the values in the image matrix. Component weights is an image matrix with the sum of the individual weights for each cell. If no weights are specified the default for each observation is one and so the weights will be the number of observations in each bin.

See Also

image.smooth, image.plot, Krig.discretize, Krig.replicates

Examples

# convert precip data to 50X50 image  
look<- as.image( RMprecip$y, x= RMprecip$x, nx=50, ny=50)
image.plot( look) 

# reduced grid extent compared to the domain
gridList<- list( x = seq(-105,-101,length.out=10),
                 y = seq(  38, 42,length.out=10) )
look2<- as.image( RMprecip$y, x= RMprecip$x,grid=gridList)
image.plot( look2) 

# number of obs in each cell -- in this case equal to the 
# aggregated weights because each obs had equal weight in the call

image.plot( look$x ,look$y, look$weights, col=terrain.colors(50)) 
# hot spot is around Denver

Description

Reformats the vector from evaluating a function on a grid of points into a list for use with surface plotting function. The list has the usual components x,y and z and is suitable for use with persp, contour, image and image.plot.

Usage

as.surface(obj, z, location=NULL, order.variables="xy")

Arguments

Details

This function was written to simply to go back and forth between a matrix of gridded values and the stacked vector obtained by stacking columns. The main application is evaluating a function at each grid point and then reforming the results for plotting. (See example below.)

If zimage is matrix of values then the input vector is c( zimage). To go from the stacked vector to the matrix one needs the the nrow ncol and explains why grid information must also be specified.

Note that the z input argument must be in the order values in order of stacking columns of the image. This is also the order of the grid points generated by make.surface.grid.

To convert irregular 2-d data to a surface object where there are missing cells see the function as.image.

Value

A list of class surface. This object is a modest generalization of the list input format (x,y,z,) for the S functions contour, image or persp.

See Also

grid.list, make.surface.grid, surface, contour, image.plot, as.image

Examples

# Make a perspective of the surface Z= X**2 -Y**2 
# Do this by evaluating quadratic function on a 25 X 25 grid
  
grid.l<-list( abcissa= seq( -2,2,,15), ordinate= seq( -2,2,,20)) 
xg<-make.surface.grid( grid.l)
# xg is a 300X2 matrix that has all pairs of X and Y grid values 
z<- xg[,1]**2 - xg[,2]**2  
# now fold z in the matrix format needed for persp 
out.p<-as.surface( xg, z) 
persp( out.p) 
# also try  plot( out.p) to see the default plot for a surface object

Description

The BD data frame has 89 rows and 5 columns. There are 89 runs with four buffer components (KCL, MgCl2, KP04, dnTP) systematically varied in a space-filliing design. The response is the DNA amplification rate.

Format

This data frame contains the following columns:

KCl

Buffer component.

MgCl2

Buffer component.

KPO4

Buffer component.

dNTP

Buffer component, deoxyribonucleotides.

lnya

Exponential amplification rate on a log scale, i.e. the actual amplification rate.

Source

Thanks to Perry Haaland and Michael OConnell.

Becton Dickinson Research Center Research Triangle Park, NC

See Also

Tps

Examples

# fitting a DNA strand 
# displacement amplification  surface to various buffer compositions 
fit<- Tps(BD[,1:4],BD$lnya,scale.type="range") 
surface(fit)  # plots fitted surface and contours

Description

Plots boxplots of several groups of data and allows for placement at different horizontal or vertical positions or colors. It is also flexible in the input object, accepting either a list or matrix.

Usage

bplot(x, by, pos=NULL, at = pos, add = FALSE, boxwex =
                 0.8,xlim=NULL, ...)

Arguments

Details

This function was created as a complement to the usual S/R function for boxplots. The current function makes it possible to put the boxplots at unequal x or y positions in a rational way using the at or pos arguments. This is useful for visually grouping a large set of boxplots into several groups. Also placement of the boxplots with respect to the axis can add information to the plot. Another aspect is the emphasis on data structures for groups of data. One useful feature is the by option to break up the x vector into distinct groups.

Use axis(3) (axis(4)) to add an axis along the top (right side) or omit the category names and draw on the bottom axis(1) (left side axis(2)).

The older bplot function drew the boxplots from scratch and if one needs to do this refer to the old functions: describe.bplot, draw.bplot.obj, bplot.xy, bplot.obj

Finally to bin data into groups based on a continuous variable and to make bplots of each group see bplot.xy.

See Also

bplot.xy

Examples

#
set.seed(123)
temp<- matrix( rnorm(12*8), ncol=12)
pos<- c(1:6,9, 12:16)*100
bplot(temp)
#
par(las=2)
bplot( temp, pos=pos, names=paste( "Data",1:12, sep=""))
# add an axis along top for reference
axis(3)

#
# Xmas boxplots in pleasing red and green 
bplot( temp, pos=pos,  col=c("red4", "green4"))
# add an axis on top
axis( 3)

Description

Draws boxplots for y by binning on x. This gives a coarse, but quick, representation of the conditional distrubtion of [Y|X] in terms of boxplots.

Usage

bplot.xy(x, y, N = 10, breaks = pretty(x, N, eps.correct = 1),
                 plot = TRUE, axes = TRUE, ...)

Arguments

See Also

bplot, draw.bplot

Examples

# condition on swim times to see how run times vary
bplot.xy( minitri$swim, minitri$run, N=5)

# bivariate normal corr= .8
set.seed( 123)
x<-rnorm( 2000)
y<- .8*x +  sqrt( 1- .8**2)*rnorm( 200)
#
bplot.xy(x,y)
#
bplot.xy( x,y, breaks=seq( -3, 3,,25) ,
                xlim =c(-4,4), ylim =c(-4,4), col="grey80", lwd=2)
points( x,y,col=3, cex=.5)

Description

This data set used be named ozone but was changed to avoid conflict with other packages. The ChicagoO3 data is a list of components, x and y. x component is longitude and latitude position of each of the 20 Chicago monitoring stations, y is the average daily ozone values over the time period 6/3/87-8/30/87. These data are used extensively for the test scripts and simple examples. The lasting scientific value is probably minimal.

Format

This data set is a list containing the following components:

lon.lat

Longitude-latitude positions of monitoring stations.

x

An approximate Cartesian set of coordinates for the locations where the units are in miles. The origin is in the center of the locations.

y

Average daily ozone values over 1987 summer.

Source

AIRS, the EPA air quality data base.

See Also

Tps, Krig

Examples

fit<- Tps(ChicagoO3$x, ChicagoO3$y) 
# fitting a surface to ozone measurements. 
surface( fit, type="I")

Description

Simulates a stationary Gaussian random field on a regular grid with unit marginal variance. Makes use of the efficient algorithm based on the FFT know as circulant embedding.

Usage

sim.rf(obj)
circulantEmbedding(obj)
circulantEmbeddingSetup(grid, M = NULL, mKrigObject = NULL,
cov.function = "stationary.cov", cov.args = NULL, delta = NULL, ...)

Arguments

Details

The functions circulantEmbedding and circulantEmbeddingSetup are more recent fields functions, more easy to read, and recommended over sim.rf. sim.rf is limited to 2D fields while circulantEmbedding can handle any number of dimensions and has some shortcuts to be efficient for the 2D case.

The simulated field has the marginal variance that is determined by the covariance function for zero distance. Within fields the exponential and matern set this equal to one ( e.g. Matern(0) ==1) so that one simulates a random field with a marginal variance of one. For stationary.cov the marginal variance is whatever Covariance(0) evaluates to and we recommend that alternative covariance functions also be normalized so that this is one.

Of course if one requires a Gaussian field with different marginal variance one can simply scale the result of this function. See the third example below.

Both sim.rf and circulantEmbedding take an object that includes some preliminary calculations and so is more efficient for simulating more than one field from the same covariance.

The algorithm using an FFT known as circulant embedding, may not always work if the correlation range is large. Specifically the weight function obtained from the FFT of the covariance field will have some negative values. A simple fix is to increase the size of the domain so that the correlation scale becomes smaller relative to the extent of the domain. Increasing the size can be computationally expensive, however, and so this method has some limitations. But when it works it is an exact simulation of the random field.

For a stationary model the covariance object ( or list) for circulantEmbedding should have minmally, the components: That is names( obj) should give "m" "grid" "M" "wght"

where m is the number of grid points in each dimension, grid is a list with components giving the grid points in each coordinate. M is the size of the larger grid that is used for "embedding" and simulation. Usually M = 2*m and results in an exact simulation of the stationary Gaussian field. The default if M is not passed is to find the smallest power of 2 greater than 2*m. wght is an array from the FFT of the covariance function with dimensions M. Keep in mind that for the final results only the array that is within the indices 1: m[i] for each dimension i is retained. This can give a much larger intermediate array, however, in the computation. E.g. if m[1] = 100 and m[2]=200 by default then M[1] = 256 and M[2] = 512. A 256 X 512 array is simluated with to get the 100 by 200 result.

The easiest way to create the object for simulation is to use circulantEmbeddingSetup.

For the older function sim.rf one uses the image based covariance functions with setup=TRUE to create the list for simulation. See the example below for this usage.

The classic reference for this algorithm is Wood, A.T.A. and Chan, G. (1994). Simulation of Stationary Gaussian Processes in [0,1]^d . Journal of Computational and Graphical Statistics, 3, 409-432. Micheal Stein and Tilman Gneiting have also made some additional contributions to the algortihms and theory.

Value

sim.rf: A matrix with the random field values.

circulantEmbedding: An array according to the grid values specified in the setup.

circulantEmbeddingetup: A list with components

"m" "grid" "dx" "M" "wght" "call"

With the information needed to simulate the field.

See Also

stationary.cov, stationary.image.cov

Examples

#Simulate a Gaussian random field with an exponential covariance function,  
#range parameter = 2.0 and the domain is  [0,5]X [0,5] evaluating the 
#field at a 100X100 grid.  
  grid<- list( x= seq( 0,5,,100), y= seq(0,5,,100)) 
  obj<- circulantEmbeddingSetup( grid, Covariance="Exponential", aRange=.5)
  set.seed( 223)
  look<-  circulantEmbedding( obj)
# Now simulate another ... 
  look2<- circulantEmbedding( obj)
# take a look at first two  
 set.panel(2,1)
 image.plot( grid[[1]], grid[[2]], look) 
 title("simulated gaussian fields")
 image.plot( grid[[1]], grid[[2]], look2) 
 title("another realization ...")
 
# Suppose one requires an exponential, range = 2
# but marginal variance = 10 ( sigma in fields notation)
look3<- sqrt( 10)*circulantEmbedding( obj)

## Not run: 
# an interesting 3D field

grid<- list(  1:40,  1:40, 1:16  )

obj<- circulantEmbeddingSetup( grid,
                         cov.args=list( Covariance="Matern", aRange=2, smoothness=1.0)
                         )
# NOTE: choice of aRange is close to giving a negative weight array
set.seed( 122)
look<- circulantEmbedding( obj )
# look at slices in the 3rd dimension 
set.panel( 4,4)
zr<- range( look)
par( mar=c(1,1,0,0))
for(  k in 1:16){
image( grid[[1]], grid[[2]], look[,,k], zlim= zr, col=tim.colors(256),
       axes=FALSE, xlab="", ylab="")
}



## End(Not run)


# same as first example using the older sim.rf

grid<- list( x= seq( 0,10,length.out=100) , y= seq( 0,10,length.out=100) )
obj<-Exp.image.cov( grid=grid, aRange=.75, setup=TRUE)
set.seed( 223)
look<- sim.rf( obj)
# Now simulate another ... 
look2<- sim.rf( obj)

Description

This is an example of moderately large spatial data set and consists of simulated CO2 concentrations that are irregularly sampled from a lon/lat grid. Also included is the complete CO2 field (CO2.true) used to generate the synthetic observations.

Usage

data(CO2)

Format

The format of CO2 is a list with two components:

• lon.lat: 26633x2 matrix of the longitude/latitude locations. These are a subset of a larger lon/lat grid (see example below).

• y: 26633 CO2 concentrations in parts per million.

The format of CO2.true is a list in "image" format with components:

• x longitude grid values.

• y latitude grid values.

• z an image matrix with CO2 concentration in parts per million

• mask a logical image that indicates with grid locations were selected for the synthetic data set CO2.

Details

This data was generously provided by Dorit Hammerling and Randy Kawa as a test example for the spatial analysis of remotely sensed (i.e. satellite) and irregular observations. The synthetic data is based on a true CO2 field simulated from a geophysical, numerical model.

Examples

## Not run: 

data(CO2)
#
# A quick look at the observations with world map
quilt.plot( CO2$lon.lat, CO2$y)
world( add=TRUE)

# Note high concentrations in Borneo (biomass burning), Amazonia and
# ... Michigan (???).

# spatial smoothing using the wendland compactly supported covariance
# see help( fastTps) for details
# First smooth using locations and Euclidean distances 
# note taper is in units of degrees 
out<-fastTps( CO2$lon.lat, CO2$y, aRange=4, lambda=2.0) 
#summary of fit note about 7300 degrees of freedom 
# associated with fitted surface
 print( out)
# image plot on a grid  (this takes a while)
surface( out, type="I", nx=300, ny=150)
# smooth with respect to great circle distance 
out2<-fastTps( CO2$lon.lat, CO2$y, lon.lat=TRUE,lambda=1.5, aRange=4*68) 
print(out2)
#surface( out2, type="I", nx=300, ny=150)

# these data are actually subsampled from a grid. 
# create the image object that holds the data
#

temp<- matrix( NA, ncol=ncol(CO2.true$z), nrow=nrow(CO2.true$z))
temp[ CO2.true$mask] <- CO2$y

# look at gridded object. 
 image.plot(CO2.true$x,CO2.true$y, temp)

# to predict _exactly_ on this grid for the second fit;
# (this takes a while)
look<- predictSurface( out2, list( x=CO2.true$x, y=CO2.true$y) )
image.plot(look)


## End(Not run)

Description

Source: These is a group of R data sets for monthly min/max temperatures and precipitation over the period 1895-1997. It is a subset extracted from the more extensive US data record. Temperature is in degrees C and precipitation is total monthly accumulation in millimeters. Note that minimum (maximum) monthly tempertuare is the mean of the daily minimum (maximum) temperatures.

Data domain:

A rectagular lon/lat region [-109.5,-101]x [36.5,41.5] larger than the boundary of Colorado comprises approximately 400 stations. Although there are additional stations reported in this domain, stations that only report preicipitation or only report temperatures have been excluded. In addition stations that have mismatches between locations and elevations from the two meta data files have also been excluded. The net result is 367 stations that have colocated temperatures and precipitation.

Format

This group of data sets is organized with the following objects:

CO.info

A data frame with columns: station id, elev, lon, lat, station name

CO.elev

elevation in meters

CO.elevGrid

An image object being elevation in meters on a 4 km grid covering Colorado.

CO.id

alphanumeric station id codes

CO.loc

locations in lon/lat

CO.Grid

Just the grid.list used in the CO.elevGrid.

CO.ppt CO.tmax CO.tmin

Monthly means as three dimensional arrays ( Year, Month, Station). Temperature is in degrees C and precipitation in total monthly accumulation in millimeters.

CO.ppt.MAM CO.tmax.MAM CO.tmin.MAM

Spring seasonal means (March, April,May) as two dimensional arrays (Year, Station).

CO.MAM.ppt.climate CO.MAM.tmax.climate CO.MAM.tmin.climate

Spring seasonal means (March, April,May) means by station for the period 1960-1990. If less than 15 years are present over this period an NA is recorded. No detreding or other adjustments have been made for these mean estimates.

Creation of data subset

Here is the precise R script used to create this data subset from the larger US monthly data set. This parent, R binary file can be obtained by contacting Doug Nychka ([email protected]).

These technical details are not needed for casual use of the data – skip down to examples for some R code that summarizes these data.

attach("RData.USmonthlyMet.bin")

#To find a subset that covers Colorado (with a bit extra):


indt<- UStinfo$lon< -101 &  UStinfo$lon > -109.5
indt<- indt & UStinfo$lat<41.5 &  UStinfo$lat>36.5

# check  US(); points( UStinfo[indt,3:4])

#find common names restricting choices to the temperature names
tn<-  match(  UStinfo$station.id, USpinfo$station.id)
indt<- !is.na(tn) & indt

# compare  metadata locations and elevations. 
# initial matches to precip stations
CO.id<- UStinfo[indt,1]
CO.names<- as.character(UStinfo[indt,5])
pn<-  match(  CO.id, USpinfo$station.id)

loc1<- cbind( UStinfo$lon[indt], UStinfo$lat[indt], UStinfo$elev[indt])
loc2<- cbind( USpinfo$lon[pn], USpinfo$lat[pn], USpinfo$elev[pn])

abs(loc1- loc2) ->  temp
indbad<- temp[,1] > .02 | temp[,2]> .02 |  temp[,3] > 100

# tolerance at 100 meters set mainly to include the CLIMAX station
# a high altitude station. 

data.frame(CO.names[ indbad], loc1[indbad,], loc2[indbad,], temp[indbad,] )


#  CO.names.indbad.      X1    X2   X3    X1.1  X2.1 X3.1 X1.2 X2.2 X3.2
#1     ALTENBERN    -108.38 39.50 1734 -108.53 39.58 2074 0.15 0.08  340
#2     CAMPO 7 S    -102.57 37.02 1311 -102.68 37.08 1312 0.11 0.06    1
#3     FLAGLER 2 NW -103.08 39.32 1519 -103.07 39.28 1525 0.01 0.04    6
#4     GATEWAY 1 SE -108.98 38.68 1391 -108.93 38.70 1495 0.05 0.02  104
#5     IDALIA       -102.27 39.77 1211 -102.28 39.70 1208 0.01 0.07    3
#6     KARVAL       -103.53 38.73 1549 -103.52 38.80 1559 0.01 0.07   10
#7     NEW RAYMER   -103.85 40.60 1458 -103.83 40.58 1510 0.02 0.02   52

# modify the indt list to exclude these mismatches (there are 7 here)

badones<-  match(  CO.id[indbad], UStinfo$station.id)
indt[ badones] <- FALSE

###### now have working set of CO stations have both temp and precip 
##### and are reasonably close to each other. 

N<- sum( indt)
# put data in time series order instead of table of year by month.
CO.tmax<-  UStmax[,,indt]
CO.tmin<-  UStmin[,,indt]

CO.id<- as.character(UStinfo[indt,1])
CO.elev<- UStinfo[indt,2]
CO.loc <-  UStinfo[indt,3:4]
CO.names<- as.character(UStinfo[indt,5])

CO.years<- 1895:1997

# now find precip stations that match temp stations
pn<-  match(  CO.id, USpinfo$station.id)
# number of orphans
sum( is.na( pn))

pn<- pn[ !is.na( pn)]
CO.ppt<- USppt[,,pn]

# checks --- all should zero

ind<- match( CO.id[45], USpinfo$station.id)
mean( abs( c(USppt[,,ind])   - c(CO.ppt[,,45]) ) , na.rm=TRUE)

ind<- match( CO.id[45], UStinfo$station.id)
mean( abs(c((UStmax[,,ind])) - c(CO.tmax[,,45])), na.rm=TRUE)

mean( abs(c((UStmin[,,ind])) - c(CO.tmin[,,45])), na.rm=TRUE)


# check order
ind<- match( CO.id, USpinfo$station.id)
 sum( CO.id != USpinfo$station.id[ind])
ind<- match( CO.id, UStinfo$station.id)
 sum( CO.id != UStinfo$station.id[ind])


# (3 4 5) (6 7 8) (9 10 11)  (12 1 2)
N<- ncol( CO.tmax)

CO.tmax.MAM<- apply( CO.tmax[,3:5,],c(1,3), "mean")

CO.tmin.MAM<- apply( CO.tmin[,3:5,],c(1,3), "mean")

CO.ppt.MAM<- apply( CO.ppt[,3:5,],c(1,3), "sum")

# Now average over 1961-1990
ind<-  CO.years>=1960 & CO.years < 1990

 temp<- stats( CO.tmax.MAM[ind,])
 CO.tmax.MAM.climate<- ifelse( temp[1,] >= 15, temp[2,], NA)

 temp<- stats( CO.tmin.MAM[ind,])
 CO.tmin.MAM.climate<- ifelse( temp[1,] >= 15, temp[2,], NA)

 CO.tmean.MAM.climate<- (CO.tmin.MAM.climate + CO.tmin.MAM.climate)/2

 temp<- stats( CO.ppt.MAM[ind,])
 CO.ppt.MAM.climate<- ifelse( temp[1,] >= 15, temp[2,], NA)


save( list=c( "CO.tmax", "CO.tmin", "CO.ppt",
              "CO.id", "CO.loc","CO.years",
              "CO.names","CO.elev", 
              "CO.tmin.MAM", "CO.tmax.MAM", "CO.ppt.MAM",
              "CO.tmin.MAM.climate", "CO.tmax.MAM.climate", 
              "CO.ppt.MAM.climate", "CO.tmean.MAM.climate"), 
               file="COmonthlyMet.rda")

Examples

data(COmonthlyMet)

#Spatial plot of 1997 Spring average daily maximum temps
 quilt.plot( CO.loc,CO.tmax.MAM[103,]  )
 US( add=TRUE)
 title( "Recorded MAM max temperatures (1997)")

# min and max temperatures against elevation

matplot( CO.elev, cbind( CO.tmax.MAM[103,], CO.tmin.MAM[103,]),
  pch="o", type="p",
  col=c("red", "blue"), xlab="Elevation (m)", ylab="Temperature (C)")
title("Recorded MAM max (red) and min (blue) temperatures 1997")

#Fitting a spatial model:
obj<- Tps(CO.loc,CO.tmax.MAM.climate,  Z= CO.elev )
## Not run: 
out<- spatialProcess(CO.loc,CO.tmax.MAM.climate, 
          smoothness=1.0, Z= CO.elev)
surface( out)
	  
## End(Not run)

Description

Adds one or more color scales in a horizontal orientation, vertical orientation to an existing plot.

Usage

colorbar.plot(x, y, strip, strip.width = 0.1, strip.length = 4 * strip.width, 
zrange = NULL, adj.x = 0.5, adj.y = 0.5, col = tim.colors(256), 
horizontal = TRUE, ...)

Arguments

Details

This function draws the strips as a sequence of image plots added to the existing plot. The main work is in creating a grid ( x,y) for the image that makes sense when superimposed on the plot. Note that although the columns of strip are considered as separate strips these can be oriented either horizontally or vertically based on the value of horizontal. The rows of zrange are essentially the zlim argument passed to the image function when each strip is drawn.

Don't forget to use locator to interactively determine positions. text can be used to label points neatly in conjunction with setting adj.x and adj.y. Although this function is inefficient for placing images at arbitrary locations on a plot the code can be easily adapted to do this.

This function was created to depict univariate posterior distribution on a map. The values are quantiles of the distribution and the strips when added under a common color scale give an overall impression of location and scale for several distributions.

Author(s)

Doug Nychka

See Also

image.plot, arrow.plot, add.image

Examples

# set up a plot but don't plot points  and no "box"
plot( 1:10, (1:10)*10, type="n", bty="n") 
# of course this could be anything 

y<- cbind( 1:15, (1:15)+25)

colorbar.plot( 2.5, 30, y)
points( 2.5,30, pch="+", cex=2, adj=.5)
# note that strip is still in 1:8 aspect even though plot has very 
# different ranges for x and y. 

# adding legend using image.plot
zr<- range( c( y))
image.plot( legend.only=TRUE, zlim= zr) 
# see help(image.plot) to create more room in margin etc. 

zr<- rbind( c(1,20), c(1,100)) # separate ranges for columns of y. 
colorbar.plot( 5, 70, y, adj.x=0, zrange= zr)
# some reference lines to show placement
xline( 5, lty=2) # strip starts at x=5 
yline(70, lty=2)  # strip is centered around y=7 (because adj.y=.5 by default)

# many strips on common scale.

y<- matrix( 1:200, ncol=10)
colorbar.plot( 2, 75, y, horizontal=FALSE, col=rainbow(256)) 

# Xmas strip
y<- cbind( rep( c(1,2),10))
y[15] <- NA # NA's should work 
colorbar.plot( 6, 45, y, adj.y=1,col=c("red", "green"))
text(6,48,"Christmas strip", cex=2)

# lennon thumbnail
# there are better ways to this ... see add.image for example.
data( lennon)
colorbar.plot( 7.5,22, lennon, 
           strip.width=.25, strip.length=.25, col=grey(seq( 0,1,,256)))

Description

compactToMat transforms a matrix from compact, vector form to a standard matrix. Only symmetric matrices can be stored in this form, since a compact matrix is stored as a vector with elements representing the upper triangle of the matrix. This function assumes the vector does not contain diagonal elements of the matrix.

An example of a matrix stored in compact form is any matrix generated from the rdist function with compact=TRUE.

Usage

compactToMat(compactMat, diagVal=0, lower.tri=FALSE, upper.tri=TRUE)

Arguments

Value

The standard form matrix represented by the input compact matrix

Author(s)

John Paige

See Also

rdist, link{dist}

Examples

################
#Calculate distance matrix from compact form:
################

#make a distance matrix
distOut = rdist(1:5, compact=TRUE)
print(distOut)

#note that distOut is in compact form:
print(c(distOut))

#convert to standard matrix form:
distMat = compactToMat(distOut)

################
#fast computation of covariance matrix:
################

#generate 5 random points on [0,1]x[0,1] square
x = matrix(runif(10), nrow=5)

#get compact distance matrix
distOut = rdist(x, compact=TRUE)

#evaluate Exponential covariance with range=1.  Note that
#Covariance function is only evaluated over upper triangle
#so time is saved.
diagVal = Exponential(0, range=1)
compactCovMat = Exponential(distOut, range=1)
upperCovMat = compactToMat(compactCovMat, diagVal)
lowerCovMat = compactToMat(compactCovMat, diagVal, lower.tri=TRUE, upper.tri=FALSE)
fullCovMat = compactToMat(compactCovMat, diagVal, lower.tri=TRUE, upper.tri=TRUE)
compactCovMat
lowerCovMat
upperCovMat
fullCovMat

Description

Given two sets of locations these functions compute the cross covariance matrix for some covariance families. In addition these functions can take advantage of spareness, implement more efficient multiplcation of the cross covariance by a vector or matrix and also return a marginal variance to be consistent with calls by the Krig function.

stationary.cov and Exp.cov have additional arguments for precomputed distance matrices and for calculating only the upper triangle and diagonal of the output covariance matrix to save time. Also, they support using the rdist function with compact=TRUE or input distance matrices in compact form, where only the upper triangle of the distance matrix is used to save time.

Note: These functions have been been renamed from the previous fields functions using 'Exp' in place of 'exp' to avoid conflict with the generic exponential function (exp(...))in R.

Usage

Exp.cov(x1, x2 = NULL, aRange = 1, p = 1, distMat = NA, C =
                 NA, marginal = FALSE, onlyUpper = FALSE, theta = NULL,
                 ...)
                 
Exp.simple.cov(x1, x2, aRange =1, C=NA,marginal=FALSE, theta=NULL)

Rad.cov(x1, x2, p = 1, m=NA, with.log = TRUE, with.constant = TRUE, 
               C=NA,marginal=FALSE, derivative=0)

cubic.cov(x1, x2, aRange = 1, C=NA, marginal=FALSE, theta=NULL) 

Rad.simple.cov(x1, x2, p=1, with.log = TRUE, with.constant = TRUE, 
               C = NA, marginal=FALSE)

stationary.cov(x1, x2=NULL, Covariance = "Exponential", Distance = "rdist", 
  Dist.args = NULL, aRange = 1, V = NULL, C = NA, marginal = FALSE, 
  derivative = 0, distMat = NA, onlyUpper = FALSE, theta=NULL, ...)

stationary.taper.cov(x1, x2, Covariance="Exponential", 
           Taper="Wendland", 
           Dist.args=NULL, Taper.args=NULL, 
           aRange=1.0,V=NULL, C=NA, marginal=FALSE,
           spam.format=TRUE,verbose=FALSE, theta=NULL,...)
	   
Tps.cov(x1, x2 = NULL, cardinalX, m=2,
                   C = NA, aRange=NA,
                   marginal = FALSE
          ) 	   

wendland.cov(x1, x2, aRange = 1, V=NULL, k = 2, C = NA, 
             marginal =FALSE,Dist.args = list(method = "euclidean"), 
             spam.format = TRUE, derivative = 0, verbose=FALSE, theta=NULL)
Paciorek.cov(x1,
                           x2 = NULL,
                           Distance = "rdist",
                           Dist.args = NULL,
                           aRangeObj = 1,
                           rhoObj = NULL,
                           C = NA,
                           marginal = FALSE,
                           smoothness = .5)

Arguments

Details

For purposes of illustration, the function Exp.cov.simple is provided in fields as a simple example and implements the R code discussed below. List this function out as a way to see the standard set of arguments that fields uses to define a covariance function. It can also serve as a template for creating new covariance functions for the Krig and mKrig functions. Also see the higher level function stationary.cov to mix and match different covariance shapes and distance functions.

A common scaling for stationary covariances: If x1 and x2 are matrices where nrow(x1)=m and nrow(x2)=n then this function will return a mXn matrix where the (i,j) element is the covariance between the locations x1[i,] and x2[j,]. The exponential covariance function is computed as exp( -(D.ij)) where D.ij is a distance between x1[i,] and x2[j,] but having first been scaled by aRange. Specifically if aRange is a matrix to represent a linear transformation of the coordinates, then let u= x1%*% t(solve( aRange)) and v= x2%*% t(solve(aRange)). Form the mXn distance matrix with elements:

D[i,j] = sqrt( sum( ( u[i,] - v[j,])**2 ) ).

and the cross covariance matrix is found by exp(-D). The tapered form (ignoring scaling parameters) is a matrix with i,j entry exp(-D[i,j])*T(D[i,j]). With T being a positive definite tapering function that is also assumed to be zero beyond 1.

Note that if aRange is a scalar then this defines an isotropic covariance function and the functional form is essentially exp(-D/aRange).

Implementation: The function r.dist is a useful FIELDS function that finds the cross Euclidean distance matrix (D defined above) for two sets of locations. Thus in compact R code we have

exp(-rdist(u, v))

Note that this function must also support two other kinds of calls:

If marginal is TRUE then just the diagonal elements are returned (in R code diag( exp(-rdist(u,u)) )).

If C is passed then the returned value is exp(-rdist(u, v)) %*% C.

Some details on particular covariance functions:

Stationary covariance stationary.cov:

Here the computation is to apply the function Covariance to the distances found by the Distance function. For example

Exp.cov(x1,x2, aRange=MyTheta)

and

stationary.cov( x1,x2, aRange=MyTheta, Distance= "rdist", Covariance="Exponential")

are the same. This also the same as

stationary.cov( x1,x2, aRange=MyTheta, Distance= "rdist", Covariance="Matern",smoothness=.5).

Radial basis functions (Rad.cov:

The functional form is Constant* rdist(u, v)**p for odd dimensions and Constant* rdist(u,v)**p * log( rdist(u,v) ) For an m th order thin plate spline in d dimensions p= 2*m-d and must be positive. The constant, depending on m and d, is coded in the fields function radbas.constant. This form is only a generalized covariance function – it is only positive definite when restricted to linear subspace. See Rad.simple.cov for a coding of the radial basis functions in R code.

Tps.cov

This covariance can be used in a standard "Kriging" computation to give a thin-plate spline (TPS). This is useful because one can use the high level function spatialProcess and supporting functions for the returned object, including conditional simulation. The standard computation for a TPS uses the radial basis functions as given in Rad.cov and uses a QR decomposition based a polynoimial matrix to reduce the dimension of the radial basis function and yield a positive definite matrix. This reduced matrix is then used in the regular compuations to find the spatial estimate. The function Krig and specifically Tps implements this algoritm. The interested reader should look at Krig.engine.default and specifically at the TMatrix polynomial matrix and resulting reduced positive definite matrix tempM. The difficulty with this approach is that is not amenable to taking advantage of sparsity in the covariance matrix.

An alternative that is suggested by Grace Wahba in Spline models for observational data is to augment the radial basis functions with a low rank set of functions based on a low order polynomial evaluated at a set of points. The set of locatios for this modifications are called cardinal points due the to property mentioned below. This is implemented in Tps.cov leading to a full rank (non-stationary!) covariance function. If the fixed part of the spatial model also includes this same order polynomial then the final result gives a TPS and is invariant to the choice of cardinal points. To streamline using this covariance when it isspecified in the spatialProcess function the cardinal points will choosen automaitcally based on the observation locations and the spline order, m using a space filling design. A simple example with fixed smoothing parameter, lambda <- .1 may help

data( "ozone2")
s<- ozone2$lon.lat
y<- ozone2$y[16,]
 
fitTps1<- spatialProcess( s,y, cov.function= "Tps.cov", lambda=.1)               
 

and compare the results to the standard algorithm.

fitTps2<- Tps( s,y, scale.type ="unscaled", lambda=.1)

stats( abs(fitTps2$fitted.values - fitTps1$fitted.values)) 

Here the default choice for the order is 2 and in two dimensions implies a linear polynomial. The arguments filled in by default are shown below

 fitTps1$args
$cardinalX
        [,1]   [,2]
[1,] -85.289 40.981
[2,] -90.160 38.330
[3,] -91.229 43.812
attr(,"scaled:scale")
[1] 1 1
attr(,"scaled:center")
[1] 0 0

$aRange
[1] NA

fitTps1$mKrig.args
[[1]]
NULL

$m
[1] 2

$collapseFixedEffect
[1] TRUE

$find.trA
[1] TRUE

cardinalX are the cardinal points chosen using a space filling criterion. These are attached to the covariance arguments list so they are used in the subsequent computations for this fit (such as predict, predictSE, sim.spatialProcess).

In general, if d is the dimension of the locations and m the order of the spline one must have 2*m-d >0. The polynomial will have choose( m+d-1,d) terms and so this many cardinal points need to be specified. As mentioned above these are chosen in a reasonable way if spatialProcess is used.

Stationary tapered covariance stationary.taper.cov :

The resulting cross covariance is the direct or Shure product of the tapering function and the covariance. In R code given location matrices, x1 and x2 and using Euclidean distance.

Covariance(rdist( x1, x2)/aRange)*Taper( rdist( x1, x2)/Taper.args$aRange)

By convention, the Taper function is assumed to be identically zero outside the interval [0,1]. Some efficiency is introduced within the function to search for pairs of locations that are nonzero with respect to the Taper. This is done by the SPAM function nearest.dist. This search may find more nonzero pairs than dimensioned internally and SPAM will try to increase the space. One can also reset the SPAM options to avoid these warnings. For spam.format TRUE the multiplication with the C argument is done with the spam sparse multiplication routines through the "overloading" of the %*% operator.

Nonstationary covariance function, Paciorek.cov

This implements the nonstationary model developed by Chris Paciorek and Mark Schervish that allows for a varying range parameter over space and also a varying marginal variance. This is still experimental and spatialProcess has not been generalized to fit the parameter surfaces. It can, however, be used to evaluate the model at fixed parameter surfaces. See the last example below.

This covariance works by specifying a object aRangeObj such that the generic call predict(aRangeObj, loc) will evaluate the aRange function at the locations loc. This object can be as simple a fit to local estimated aRange parameters using a fields function such as Tps or spatialProcess. More specific applications one can create a special predict function. For example suppose log aRange follows a linear model in the spatial coordinates and these coefficients are the vector Afit. Define a class and a predict function.

 

aRangeObj<- list(coef=Afit)
class(aRangeObj)<- "myclass"

predict.myclass<- function( aRangeObj, x){
aRange<- exp(cbind( 1,x) %*% aRangeObj$coef)
return( aRange)
}

Now use spatialProcess with this object and make sure predict.myclass is also loaded.

A similar strategy will also work for a varying marginal variance by creating sigmaObj and if needed a companion predict method.

About the FORTRAN: The actual function Exp.cov and Rad.cov call FORTRAN to make the evaluation more efficient this is especially important when the C argument is supplied. So unfortunately the actual production code in Exp.cov is not as crisp as the R code sketched above. See Rad.simple.cov for a R coding of the radial basis functions.

Value

If the argument C is NULL the cross covariance matrix is returned. In general if nrow(x1)=m and nrow(x2)=n then the returned matrix will be mXn. Moreover, if x1 is equal to x2 then this is the covariance matrix for this set of locations.

If C is a vector of length n, then returned value is the multiplication of the cross covariance matrix with this vector.

See Also

Krig, rdist, rdist.earth, gauss.cov, Exp.image.cov, Exponential, Matern, Wendland.cov, mKrig

Examples

# exponential covariance matrix ( marginal variance =1) for the ozone
#locations 
out<- Exp.cov( ChicagoO3$x, aRange=100)

# out is a 20X20 matrix

out2<- Exp.cov( ChicagoO3$x[6:20,],ChicagoO3$x[1:2,], aRange=100)
# out2 is 15X2 matrix 

# Kriging fit where the nugget variance is found by GCV 
# Matern covariance shape with range of 100.
# 

fit<- Krig( ChicagoO3$x, ChicagoO3$y, Covariance="Matern", aRange=100,smoothness=2)

data( ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
# Omit the NAs
good<- !is.na( y)
x<- x[good,]
y<- y[good]


# example of calling the taper version directly 
# Note that default covariance is exponential and default taper is 
# Wendland (k=2).

stationary.taper.cov( x[1:3,],x[1:10,] , aRange=1.5, Taper.args= list(k=2,aRange=2.0,
                       dimension=2) )-> temp
# temp is now a tapered 3X10 cross covariance matrix in sparse format. 

 is.spam( temp)  # evaluates to TRUE

# should be identical to
# the direct matrix product

 temp2<- Exp.cov( x[1:3,],x[1:10,], aRange=1.5) * Wendland(rdist(x[1:3,],x[1:10,]), 
                      aRange= 2.0, k=2, dimension=2)
 test.for.zero(  as.matrix(temp), temp2)

# Testing that the "V" option works as advertized ...
x1<- x[1:20,]
x2<- x[1:10,]

V<- matrix( c(2,1,0,4), 2,2)
Vi<- solve( V)

u1<- t(Vi%*% t(x1))
u2<- t(Vi%*% t(x2))

look<- exp(-1*rdist(u1,u2))
look2<- stationary.cov( x1,x2, V= V)
test.for.zero( look, look2)


# Here is an example of how the cross covariance multiply works
# and lots of options on the arguments


 Ctest<- rnorm(10)
 
 temp<- stationary.cov( x,x[1:10,], C= Ctest, 
        Covariance= "Wendland", 
            k=2, dimension=2, aRange=1.5 )

# do multiply explicitly

 temp2<- stationary.cov( x,x[1:10,],
        Covariance= "Wendland",
            k=2, dimension=2, aRange=1.5 )%*% Ctest

 test.for.zero( temp, temp2)


# use the tapered stationary version 
# cov.args is part of the argument list passed to stationary.taper.cov
# within Krig. 
# This example needs the spam package.
# 

## Not run: 

Krig(x,y, cov.function = "stationary.taper.cov", aRange=1.5,
      cov.args= list(Taper.args= list(k=2, dimension=2,aRange=2.0) )
           ) -> out2 
# NOTE: Wendland is the default taper here. 

## End(Not run)

# BTW  this is very similar to 
## Not run: 
 Krig(x,y, aRange= 1.5)-> out

## End(Not run)

##################################################
#### nonstationary covariance Paciorek.cov
##################################################
## Not run: 
M<- 20
gridList<- list(x=seq( 0,1,length.out=M),
                y=seq( 0,1,length.out=M))
sGrid<- make.surface.grid(gridList )
# An aRange surface
aRangeObj<- list(coef= c( 1,4,0))
class(aRangeObj)<- "myclass"

predict.myclass<- function( aRangeObj, x){
aRange<- exp(cbind( 1,x) %*% aRangeObj$coef)
return( aRange)
}

covMatrix<- Paciorek.cov( sGrid, sGrid, aRangeObj=aRangeObj)
# examine correlation surface between selected locations and the  full grid. 
set.panel( 2,2)
{
imagePlot( as.surface(sGrid, covMatrix[,10]))
imagePlot( as.surface(sGrid, covMatrix[,205]))
imagePlot( as.surface(sGrid, covMatrix[,305]))
imagePlot( as.surface(sGrid, covMatrix[,390]))

}

# simulation of the field
set.seed(222)
n<- nrow( sGrid)
f<- t(chol(covMatrix)) %*% rnorm(M^2)
set.panel()
imagePlot( as.surface(sGrid,f))
y<- f + .05*rnorm(n)
fitP<- spatialProcess( sGrid, y, cov.function="Paciorek.cov",
           cov.args= list(aRangeObj = aRangeObj ) )
# check estimated tau and sigma 
fitP$summary

# fitted surface
surface( fitP)


## End(Not run)

Description

Evaluates the covariance over the upper triangle of a distance matrix rather than over the entire matrix to reduce computation time. Note that the chol function only requires the upper triangle of the covariance matrix to perform the Cholesky decomposition.

Usage

ExponentialUpper(distMat, range = 1, alpha = 1/range, theta = NULL)

Arguments

Value

The covariance matrix, where only the upper triangle is calculated.

Author(s)

John Paige

See Also

Exponential

Examples

set.seed(123)

#a  distance matrix 
coords = matrix(runif(10), ncol=2)
distMat = rdist(coords)

#compute covariance matrix, but only over the upper triangle
upperCov = ExponentialUpper(distMat, range=.1)

print(distMat)
print(upperCov)

Description

Finds the set of points on a discrete grid (Candidate Set) which minimize a geometric space-filling criterion. The strength of this method is that the candidate set can satisfy whatever constraints are important for the problem.

Usage

cover.design(R, nd, nruns = 1, nn = TRUE, num.nn = 100, fixed = NULL, 
    scale.type = "unscaled", R.center, R.scale, P = -20, Q = 20,
    start = NULL, DIST = NULL, return.grid = TRUE, return.transform = 
TRUE, max.loop=20, verbose=FALSE)

Arguments

Details

OTHER DISTANCE FUNCTIONS: You can supply an R/S-function to be used as the distance metric. The expected calling sequence for this distance function is function( X1,X2){....} where X1 and X2 are matrices with coordinates as the rows. The returned value of this function should be the pairwise distance matrix. If nrow( X1)=m and nrow( X2)=n then the function should return an m by n matrix of all distances between these two sets of points. See the example for Manhattan distance below.

The candidate set and DIST function can be flexible and the last example below using sample correlation matrices is an example.

COVERAGE CRITERION: For nd design points in the set D and nc candidate points ci in the set C, the coverage criteria is defined as:

M(D,C) = [sum(ci in C) [sum(di in D) (dist(di,ci)**P]**(Q/P)]**(1/Q)

Where P, less than 0, and Q, greater than 0, are parameters. The algorithm used in "cover.design" to find the set of nd points in C that minimize this criterion is an iterative swapping algorithm which will be described briefly. The resulting design is referred to as a "coverage design" from among the class of space-filling designs. If fixed points are specified they are simply fixed in the design set and are not allowed to be swapped out.

ALGORITHM: An initial set of nd points is chosen randomly if no starting configuration is provided. The nc x nd distance matrix between the points in C and the points in D is computed, and raised to the power P. The "row sums" of this matrix are computed. Denote these as rs.i and the vector of row sums as rs. Using rs, M(D,C) is computed as:

[sum i (rs.i)**(Q/P)]**(1/Q)

Note that if point d.i is "swapped" for point c.j, one must only recompute 1 column of the original distance matrix, and 1 row. The row elements not in the ith column will be the same for all j and so only need computing when the first swapping occurs for each d.i . Denote the sum of these off-i elements as "newrow(i)". The index is i here since this is the same for all rows (j=1,...nc). Thus, for each swap, the row sums vector is updated as

rs(new) = rs(old) - column(i,old) + column(i,new)

And the jth element of rs(new) is replaced by:

rs(new)[j] = column(i,new)[k] + newrow(i)

Finally, M(D,C) is computed for this swap of the ith design point for the jth candidate point using [2]. The point in C that when swapped produces the minimum value of M(D,C) replaces d.i. This is done for all nd points in the design, and is iterated until M(D,C) does not change. When the nearest neighbor option is selected, then the points considered for swapping are limited to the num.nn nearest neighbors of the current design point.

STABILITY

The algorithm described above is guaranteed to converge. However, upon convergence, the solution is sensitive to the initial configuration of points. Thus, it is recommended that multiple optimizations be done (i.e. set nruns greater than 1 ). Also, the quality of the solution depends on the density of the points on the region. At the same time, for large regions , optimization can be computationally prohibitive unless the nearest neighbor option is employed.

Value

Returns a design object of class spatialDesign. Subscripting this object has the same effect as subscripting the first component (the design). The returned list has the following components:

References

Johnson, M.E., Moore, L.M., and Ylvisaker, D. (1990). Minimax and maximin distance designs. Journal of Statistical Planning and Inference 26, 131-148. SAS/QC Software. Volume 2: Usage and Reference. Version 6. First Edition (1995). "Proc Optex". SAS Institute Inc. SAS Campus Drive,

See Also

rdist, rdist.earth

Examples

##
## 
# first generate candidate set
set.seed(123) # setting seed so that you get the same thing I do!
test.df <- matrix( runif( 600), ncol=3)

test1.des<-cover.design(R=test.df,nd=10)

summary( test1.des)
plot( test1.des)

#

## Not run: 
candidates<- make.surface.grid( list( seq( 0,5,,20), seq(0,5,,20)))
out<- cover.design( candidates, 15)

# find 10 more points keeping this original design fixed

out3<-cover.design( candidates, 10,fixed=out$best.id)
# see what happened

plot( candidates[,1:2], pch=".")
points( out$design, pch="x")
points( out3$design, pch="o")    

# here is a strange graph illustrating the swapping history for the
# the first design. Arrows show the swap done  
# at each pass through the design.

h<- out$history
cd<- candidates
plot( cd[,1:2], pch=".")
points( out$design, pch="O", col=2)
points( out$start.design, pch="x", col=5)  

arrows(
cd[h[,2],1],
cd[h[,2],2],
cd[h[,3],1],
cd[h[,3],2],length=.1)
text( cd[h[,2],1],
cd[h[,2],2], h[,1], cex=1.0 )
                               

#
# try this out using "Manhattan distance"
#  ( distance following a grid of city streets)

dist.man<- function(x1,x2) {
            d<- ncol( x1)
            temp<- abs(outer( x1[,1], x2[,1],'-'))
            for ( k in 2:d){
               temp<- temp+abs(outer( x1[,k], x2[,k],'-'))
            }
            temp }

# use the design from the Euclidean distance as the starting
#configuration.

cover.design( candidates, 15, DIST=dist.man, start= out3$best.id)-> out2
# this takes a while ...
plot( out2$design)
points( out3$design, col=2)

# find a design on the sphere
#

candidates<- make.surface.grid( list( x=seq( -180,180,,20), y= seq( -85,
85,,20)))

out4<-cover.design( candidates, 15, DIST=rdist.earth)
# this takes a while 

plot( candidates, pch="+", cex=2)
points(out4$design, pch="o", cex=2, col="blue")

# covering based on correlation for 153 ozone stations
#
data( ozone2)

cor.mat<-cor( ozone2$y, use="pairwise")

cor.dist<- function( x1,x2)
{matrix( 1-cor.mat[ x1,x2], ncol=length(x2))}

#
# find 25 points out of the 153
# here the "locations" are just the index but the distance is 
# determined by the correlation function. 
#
out5<-cover.design(cbind(1:153),25, DIST= cor.dist, scale.type="unscaled") 

plot( ozone2$lon.lat, pch=".")
points(  ozone2$lon.lat[out5$best.id,],pch="O", col=4)
#
# this seems a bit strange probably due some funny correlation values
#

# reset panel
set.panel(1,1)

## End(Not run)

Description

Function to produce the usual wireframe perspective plot with the facets being filled with different colors. By default the colors are assigned from a color bar based on the z values. drape.color can be used to create a color matrix different from the z matrix used for the wireframe.

Usage

drape.plot(x, y, z, z2=NULL, col = tim.colors(64), zlim = range(z, na.rm=TRUE), 
 zlim2 = NULL, add.legend = TRUE, horizontal = TRUE, theta = 30, phi = 20, 
   breaks=NA, ...)

drape.color(z, col = tim.colors(64), zlim = NULL,breaks,
   transparent.color = "white", midpoint=TRUE, eps=1e-8)

Arguments

Details

The legend strip may obscure part of the plot. If so, add this as another step using image.plot.

When using drape.color just drop the results into the col argument of persp. Given this function there are no surprises how the higher level drape.plot works: it calls drape.color followed by persp and finally the legend strip is added with image.plot.

The color scales essentially default to the ranges of the z values. However, by specifying zlim and/or zlim2 one has more control of how the perspective plot is scaled and the limits of the color scale used to fill the facets. The color assignments are done by dividing up the zlim2 interval into equally spaced bins and adding a very small inflation to these limits. The mean z2 values, comprising an (M-1)X(N-1) matrix, for each facet are discretized to the bins. The bin numbers then become the indices used for the color scale. If zlim2 is not specified it is the range of the z2 matrix is used to generate the ranges of the color bar. Note that this may be different than the range of the mean facets. If z2 is not passed then z is used in its place and in this case the zlim2 or zlim argument can used to define the color scale.

This kind of plot is also supported through the wireframe function in the lattice package. The advantage of the fields version is that it uses the standard R graphics functions – and is written in R code.

The drape plot is also drawn by the fields surface function with type="P".

Value

drape.plot If an assignment is made the projection matrix from persp is returned. This information can be used to add additional 3-d features to the plot. See the persp help file for an example how to add additional points and lines using the trans3d function and also the example below.

drape.color If dim( z) = M,N this function returns a list with components:

Author(s)

D. Nychka

See Also

image.plot, quilt.plot, persp, plot.surface, surface, lattice, trans3d

Examples

# an obvious choice:
# Dr. R's favorite New Zealand Volcano!
data( volcano)
M<- nrow( volcano)
N<- ncol( volcano)
x<- seq( 0,1,,M)
y<- seq( 0,1,,N)

pm<- drape.plot( x,y,volcano, col=terrain.colors(128)) 

# use different range for color scale and persp plot
# setting of border omits the mesh lines

 drape.plot( x,y,volcano, col=topo.colors(128),zlim=c(0,300),
                     zlim2=c( 120,200), border=NA)

# note tranparent color for facets outside the zlim2 range

#The projection has been saved in pm
# add a point marking the summit
zsummit <- max( volcano) 
ix<- row( volcano)[volcano==zsummit]
iy <- col( volcano)[volcano==zsummit]
uv <- trans3d( x[ix], y[iy],zsummit,pm)
points( uv, col="magenta", pch="+", cex=2)

# overlay volcano wireframe with gradient in x direction. 

dz<- ( 
     volcano[1:(M-1), 1:(N-1)] - volcano[2:(M), 1:(N-1)] +
     volcano[1:(M-1), 2:(N)] - volcano[2:(M), 2:(N)]  
         )/2

# convert dz to a color scale:
  zlim<- range( c( dz), na.rm=TRUE)
  zcol<-drape.color( dz, zlim =zlim, col = viridis(64) )$color.index

# with these colors 

  persp( volcano, col=zcol, theta=30, phi=20,
  border=NA,expand=.3 )

# add legend using image.plot function 
  image.plot( zlim=zlim, legend.only =TRUE, horizontal =TRUE,
            col= viridis(64))

Description

This function shades the region vertically between two functions, specified as pairs of x and y vectors, and draws the functions in a darker shade. More formally, it shades all points (x,y) such that f1(x) < y < f2(x) or f2(x) < y < f1(x). When both functions have the same group of x values, the x values only need to be set once but y2 needs to be passed in by name. If the two functions intersect, the vertical space between the functions will be shaded on both sides, as implied in the definition above.

Usage

envelopePlot(x1, y1, x2 = x1, y2,
                         col ="thistle1" , lineCol = "thistle3", ...)

Arguments

Author(s)

Matt Iverson

Examples

x <- seq(0, 2*pi,, 100)
y1 <- cos(x)
y2 <- sin(x)
plot(x, y1, type="l")
envelopePlot(x, y1, y2=y2)

x1 <- c(0, 0.5, 1)
y1 <- c(0, 2, 1)
x2 <- c(0, 1)
y2 <- c(-1, 0)
plot(x1, y1, type="l", ylim=c(-1, 2))
envelopePlot(x1, y1, x2, y2)

Description

Functional form of covariance function assuming the argument is a distance between locations. As they are defined here, they are in fact correlation functions. To set the marginal variance (sill) parameter, use the sigma argument in mKrig or Krig. To set the nugget variance, use the tau2 argument in mKrig or Krig.

Usage

Exponential(d, aRange = 1, phi = 1, theta = NULL, range = NULL)
Matern(d , range = 1,alpha=1/range, smoothness = 0.5, 
       nu= smoothness, phi=1.0) 
Matern.cor.to.range(d, nu, cor.target=.5, guess=NULL,...)
RadialBasis(d,M,dimension, derivative = 0)

Arguments

Details

Exponential:

exp( -d/aRange)

Matern:

con*(d**nu) * besselK(d , nu )

Matern covariance function transcribed from Stein's book page 31 nu==smoothness, alpha == 1/range

GeoR parameters map to kappa==smoothness and phi == range check for negative distances

con is a constant that normalizes the expression to be 1.0 when d=0.

Matern.cor.to.range: This function is useful to find Matern covariance parameters that are comparable for different smoothness parameters. Given a distance d, smoothness nu, target correlation cor.target and range aRange, this function determines numerically the value of aRange so that

Matern( d, range=aRange, nu=nu) == cor.target

See the example for how this might be used.

Radial basis functions:

   C.m,d  r**(2m-d)        d- odd

   C.m,d  r**(2m-d)ln(r)    d-even

where C.m.d is a constant based on spline theory and r is the radial distance between points. See radbas.constant for the computation of the constant.

Value

For the covariance functions: a vector of covariances.

For Matern.cor.to.range: the value of the range parameter.

Author(s)

Doug Nychka

References

Stein, M.L. (1999) Statistical Interpolation of Spatial Data: Some Theory for Kriging. Springer, New York.

See Also

stationary.cov, stationary.image.cov, Wendland,stationary.taper.cov rad.cov

Examples

# a Matern correlation function 
 d<- seq( 0,10,,200)
 y<- Matern( d, range=1.5, smoothness=1.0)
 plot( d,y, type="l")

# Several Materns of different smoothness with a similar correlation 
# range

# find ranges for nu = .5, 1.0 and 2.0 
# where the correlation drops to .1 at a distance of 10 units.

 r1<- Matern.cor.to.range( 10, nu=.5, cor.target=.1)
 r2<- Matern.cor.to.range( 10, nu=1.0, cor.target=.1)
 r3<- Matern.cor.to.range( 10, nu=2.0, cor.target=.1)

# note that these equivalent ranges
# with respect to this correlation length are quite different
# due the different smoothness parameters. 

 d<- seq( 0, 15,,200)
 y<- cbind(  Matern( d, range=r1, nu=.5),
             Matern( d, range=r2, nu=1.0),
             Matern( d, range=r3, nu=2.0))

 matplot( d, y, type="l", lty=1, lwd=2)
 xline( 10)
 yline( .1)

Description

fields is a collection of functions for curve and function fitting with an emphasis on spatial data and spatial statistics. It was developed over 20+ years to provide easy to use but sophisticated tools for analyzing spatial data, particularly that encountered in the environmental sciences. For the impatient users, jump to the examples below to see how easy this is use. Please send bugs and questions to Doug Nychka, [email protected]. Positive comments are also welcome!

The major methods implemented include cubic and thin plate splines, universal Kriging and Kriging for large data sets. A more modern framework for Kriging is spatial process estimation with covariance parameters determined by maximum likelihood and the uncertainty derived from assumptions of a Gaussian process. Throughout we try to include reasonable defaults in functions that reflect our experience with analyzing spatial data. For example, the Matern covariance function is the default choice for the main spatial method.

A key feature of this package is any covariance function implemented in R code can be used for spatial prediction through the spatial functions. Another important feature is that fields will take advantage of compactly supported covariance functions in a seamless way through the spam package. See library( help=fields) for a listing of all the fields contents. We also recommend the thoughtful vignette created by Ashton Weins, Mitchell Krock, and Emma Lilly (fieldsVignette.pdf) in the fields github repository.

fields strives to have readable and tutorial code. Take a look at the source code for mKrig to see how things work "under the hood". E.g. how a linear algebra computation is overloaded to handle sparse matrices and how an output object is built up sequentially throughout a computation.

The fields source code is liberally commented. Unfortunately on loading this package, R will strip comments from the source for efficiency. You can go to CRAN fields page to download the latest "tarball" ( aka Package Source) and unzip to get code with comments. We also keep the most recent version of this package at the fields github repository. and for commented source go to the the subdirectory fields/R. This may be a more recent version, however, than what is posted to CRAN.

Details

Major methods

• spatialProcess An easy to use method that fits a spatial process model ( e.g. Kriging) but also estimates the key spatial parameters: nugget variance, sill variance and range parameter by maximum likelihood. The default covariance model is a Matern covariance. This function and related functions called by this are the core methods in fields and have much flexibility.

  spatialProcess allows one to supply a covariance function that is written in native R code. See (stationary.cov) that includes several families of covariances including the Matern and several distance metrics including great circle distance . sim.spatialProcess and simLocal.spatialProcess provide "one liners" for conditional simulation of the fitted surface.

• Tps Thin Plate spline regression including GCV and REML estimates for the smoothing parameter. For moderate size data sets as a first look we use Tps all the time. See also fastTps for an approximate method to handle very large numbers of spatial locations. Also see the help file for spatialProcess to see how to fit a thin plate plate using the more extensive set of spatial stats functions.

• sreg , splint Fast 1-D cubic smoothing splines and interpolating splines, a workhorse algorithm for more EDA and more complicated methods.

• mKrig (micro Krig) Efficient Universal Kriging and Gaussian process function, that can take advantage of sparse covariance functions and is the core algorithm called by optimization functions and for spatial predictio.

• QTps A easy to use extension of thin plate splines for quantile and robust surface fitting.

• mKrigMLEGrid and mKrigMLEJoint for maximum likelihood estimates of covariance parameters. These functions also handle replicate fields, assumed to be independent realizations, at the same locations and can also take any covariate function function written in R following the fields format

Other noteworthy functions

• vgram and vgram.matrix find variograms for spatial data (and with temporal replications.

• cover.design Generates space-filling designs where the distance function is expresed in R code.

• There are many convenient functions for working with image data and rationally (well, maybe reasonably) creating and placing a color scale on plots. This suite of tools are for the users who want to extend the "base R grahics" and retain control over details. See the last amusing example in help(imagePlot) for an example. as.image, imagePlot, bubblePlot, drape.plot, quilt.plot add.image, crop.image, half.image, average.image, designer.colors, color.scale, in.poly See also grid.list for how fields works with grids and US and world for adding a map quickly.

Generic functions that support the methods

plot - diagnostic plots of fit
summary- statistical summary of fit
print- shorter version of summary
surface- graphical display of fitted surface
predict- evaluation fit at arbitrary points
predictSE- prediction standard errors at arbitrary points.
sim.rf- Simulate a random fields on a 2-d grid.

Getting Started

Try some of the examples from help files for spatialProcess or Tps.

Some fields datasets

• CO2 Global satelite CO2 concentrations (simulated field)

• COmonthlyMet Monthly mean temperatures and precip for Colorado

• glacier An elevation dataset of a glacier also used by the applied math community to test interpolation methods.

• lennon Image of John Lennon

• NorthAmericanRainfall 50+ year average and trend for summer rainfall at 1700+ stations.

• ozone2 Daily max 8 hour ozone concentrations for the US midwest for summer 1987.

• PRISMelevation Digital elevations for the continental US at approximately 4km resolution

• rat.diet Small paired study on rat food intake over time.

• RCMexample Regional climate model output

• RMelevation Digital elevations for the Rocky Mountain Empire

• WorldBankCO2 Demographic and carbon emission data for 75 countries and for 1999.

DISCLAIMER:

The authors can not guarantee the correctness of any function or program in this package.

Examples

# some air quality data, daily surface ozone measurements for 
# the Midwest:

data(ozone2)
s<-ozone2$lon.lat
y<- ozone2$y[16,] # June 18, 1987 

# quick plot of spatial data with map
bubblePlot( s,y)
US( add=TRUE) # add US map

# fitting a thin plate spline surface (always a good place to 
# start). Here the  default smoothing (aka lambda) found by cross-validation
  fit0<- Tps(s,y)
# fits a GCV thin plate smoothing spline surface to ozone measurements.
# Hey, it does not get any easier than this!

  summary(fit0) #diagnostic summary of the fit 
  set.panel(2,2)
  plot(fit0) # four diagnostic plots of fit and residuals.

# quick plot of predicted surface
  set.panel()
  surface(fit0) # contour/image plot of the fitted surface
# see also predictSurface for more control over the evaluation grid
#
  US( add=TRUE, col="magenta", lwd=2) # US map overlaid
  title("Daily max 8 hour ozone in PPB,  June 18th, 1987")

####
  fit2<- spatialProcess( s,y)
# a "Kriging" model. The covariance defaults to a Matern 
# with smoothness 1.0.
# the nugget, sill and range parameters are found by maximum likelihood
# summary, plot, and surface also work for \code{fit2} !

  surface(fit2) # contour/image plot of the fitted surface
  US( add=TRUE, col="magenta", lwd=2) # US map overlaid
  title("Daily max 8 hour ozone in PPB,  June 18th, 1987")
## Not run: 
# And 20 approximate conditional draws of the spatial field on a grid
# with uncertainty in the 120PPB contour 
   look<- simLocal.spatialProcess(fit2, M=20)
for( k in 1:20){
contour( look$x, look$y, look$z[,,k], add=TRUE, level=c(120),
  col="white",  drawlabels=FALSE)
}


## End(Not run)

Description

Some of the basic methods in fields can be tested by directly implementing the linear algebra using matrix expressions and other functions can be cross checked within fields. These comparisons are done in the the R source code test files in the tests subdirectory of fields. The function test.for.zero is useful for comparing the tests in a meaninful and documented way.

Usage

test.for.zero( xtest, xtrue,  tol= 1.0e-8, relative=TRUE, tag=NULL)

Arguments

Details

IMPORTANT: If the R object test.for.zero.flag exists with any value ( e.g. test.for.zero.flag <- 1 ) then when the test fails this function will also generate an error in addition to printing a message. This option is added to insure that any test scripts will generate an error when any individual test fails.

An example:

> test.for.zero( 1:10, 1:10 + 1e-10, tag="First test")
Testing:  First test
PASSED test at tolerance  1e-08

> test.for.zero( 1:10, 1:10 + 1e-10, tag="First test", tol=1e-12)
Testing:  First test
FAILED test value =  1.818182e-10  at tolerance  1e-12

> test.for.zero.flag <- 1
Testing:  First test
FAILED test value =  1.818182e-10  at tolerance  1e-12
Error in test.for.zero(1:10, 1:10 + 1e-10, tag = "First test", tol = 1e-12) :

The scripts in the tests subdirectory are

Krig.test.R:

Tests basic parts of the Krig and Tps functions including replicated and weighted observations.

Krig.se.test.R:

Tests computations of standard errors for the Kriging estimate.

Krig.se.grid.test.R

Tests approximate standard errors for the Krig function found by Monte Carlo conditional simulation.

Krig.test.W.R

Tests predictions and A matrix when an off diagonal observation weight matrix is used.

Krig.se.W.R

Tests standard errors when an off diagonal observation weight matrix is used.

spam.test.R

Tests sparse matrix formats and linear algebra.

Wend.test.R

Tests form for Wendland covariance family and its use of sparse matrix formats.

diag.multiply.test.R

Tests special (efficient) version of matrix multiply for diagonal matrices.

evlpoly.test.R

Tests evaluation of univariate and multivariate polynomial evaluation.

mKrig.test.R

Tests the micro Krig function with and without sparse matrix methods.

To run the tests just attach the fields library and source the testing file. In the fields source code these are in a subdirectory "tests". Compare the output to the "XXX.Rout.save" text file.

test.for.zero is used to print out the result for each individual comparison. Failed tests are potentially bad and are reported with a string beginning

"FAILED test value = ... "

If the object test.for.zero.flag exists then an error is also generated when the test fails.

FORM OF COMPARISON: The actual test done is the sum of absolute differnces:

test value = sum( abs(c(xtest) - c( xtrue) ) ) /denom

Where denom is either mean( abs(c(xtrue))) for relative error or 1.0 otherwise.

Note the use of "c" here to stack any structure in xtest and xtrue into a vector.

Description

Some supporting functions that are internal to fields top level methods. Variants of these might be found in the R base but these have been written for cleaner code or efficiency.

Usage

fields.duplicated.matrix(mat, digits = 8) 

fields.mkpoly(x, m = 2, tag = "term")

fields.derivative.poly(x, m,dcoef)

fields.evlpoly( x, coef)

fields.evlpoly2( x, coef, ptab)

Arguments

Details

fields.duplicated finds duplicate rows in a matrix. The digits arguments is the number of digits that are considered in the comparison. The returned value is an array of integers from 1:M where M is the number of unique rows and duplicate rows are referenced in the same order that they appear as the rows of mat.

fields.mkpoly computes the complete matrix of all monomial terms up to degree (m-1). Each row of x is are the componets of a vector. (The fields function mkpoly returns the number of these terms.) In 2 dimensions with m=3 there 6 polynomial terms up to quadratic ( 3-1 =2) order and will be returned as the matrix:

cbind( 1 , x[,1], x[,2], x[,1]**2, x[,1]*x[,2], x[,2]**2 )

This function is used for the fixed effects polynomial or spatial drift used in spatial estimating functions Krig, Tps and mKrig. The matrix ptab is a table of the powers in each term for each variable and is included as an attribute to the matrix returned by this function. See the attr function for extracting an attribute from an object.

ptab for the example above is

    [,1] [,2]
[1,]    0    0
[2,]    1    0
[3,]    0    1
[4,]    2    0
[5,]    1    1
[6,]    0    2

This information is used in finding derivatives of the polynomial is also used to create column names for the terms that are of higher order than linear.

fields.deriviative.poly finds the partial derivative matrix of a multidimensional polynomial of degree (m-1) at different vector values and with coefficients dcoef. This function has been orgainzed to be a clean utility for the predicting the derivative of the estimated function from Krig or mKrig. Within the fields context the polynomial itself would be evaluated as fields.mkpoly( x,m)%*%dcoef. If x has d columns ( also the dimension of the polynomial) and n rows the partial derivatives of this polynomial at the locations x can be organized in a nXd matrix. This is the object returned by ths function.

evlpoly and evlpoly2 are FORTRAN based functions for evaluating univariate polynomials and multivariate polynomials. The table of powers (ptab) needed for evlpoly2 is the same format as that returned my the fields.mkpoly function.

Author(s)

Doug Nychka

Description

This is an extended example for using the sparse/fast interpolation methods in mKrig to evaluate a Kriging estimate on a large grid.

Details

mKrig is a flexible function for surface fitting using a spatial process model. It can also exploit sparse matrix methods forlarge data sets by using a compactly supported covariance. The example below shows how ot evaluate a solution on a big grid. (Thanks to Jan Klennin for this example.)

Examples

x<- RMprecip$x
y<- RMprecip$y

Tps( x,y)-> obj

# make up an 80X80 grid that has ranges of observations
# use same coordinate names as the x matrix

glist<- fields.x.to.grid(x, nx=80, ny=80) # this is a cute way to get a default grid that covers x

# convert grid list to actual x and y values ( try plot( Bigx, pch="."))
    make.surface.grid(glist)-> Bigx 

# include actual x locations along with grid. 
    Bigx<- rbind( x, Bigx)

# evaluate the surface on this set of points (exactly)

    predict(obj, x= Bigx)-> Bigy

# set the range for the compact covariance function 
# this will involve  less than 20 nearest neighbors that have
# nonzero covariance
# 

     V<- diag(c( 2.5*(glist$lon[2]-glist$lon[1]), 
                 2.5*(glist$lat[2]-glist$lat[1])))
## Not run: 
# this is an interplotation of the values using a compact 
# but thin plate spline like covariance. 
    mKrig( Bigx,Bigy, cov.function="wendland.cov",k=4, V=V, 
                 lambda=0)->out2 
# the big evaluation this takes about 45 seconds on a Mac G4 latop
    predictSurface( out2, nx=400, ny=400)-> look

## End(Not run)

# the nice surface
## Not run:     
    surface( look)
    US( add=TRUE, col="white")

## End(Not run)

Description

Here are some technical hints for assembling multiple plots with common legends or axes and setting the graphics parameters to make more readable figures. Also we an index to the defaultcolors in R graphics and setting their definitions in LaTeX. All these hints use the standard graphics environment.

Usage

fields.style()
fields.color.picker()

Details

fields.style is a simple function to enlarge the characters in a plot and set the colors. List this out to modify the choices.

##Two examples of concentrating a panel of plots together
## to conserve the white space. 
## see also the example in image.plot using split.screen. 
## The basic trick is to use the oma option to reserve some space around the 
## plots.  Then unset the outer margins to use that room. 

library( fields)

# some hokey image data
x<- 1:20
y<- 1:15
z<-  outer( x,y,"+")
zr<- range( c(z))

# add common legend to 3X2 panel 

par( oma=c(4,0,0,0))
set.panel( 3,2)
par( mar=c(1,1,0,0))

# squish plots together with just 1 space between
for ( k in 1:6){
image( x,y,z, axes=FALSE, xlab="", ylab="", zlim=zr)
}

par( oma=c(0,0,0,0))
image.plot( zlim=zr, legend.only=TRUE, horizontal=TRUE, legend.mar=5)

# you may have to play around with legend.mar and the oma settings to
# get enough space.



##
### also add some axes on the sides. in a lattice style 
## note oma adds some more room at bottom. 

par( oma=c(8,6,1,1))
set.panel( 3,2)
par( mar=c(1,1,0,0))
##
for ( k in 1:6){
 image( x,y,z, axes=FALSE, xlab="", ylab="", zlim=zr)
 box() # box around figure

# maybe draw an x axis 
  if( k  %in% c(5,6) ){
  axis( 1, cex.axis=1.5)
  mtext( line=4, side=1, "Xstuff")}

# maybe draw a y axis
  if( k  %in% c(1,3,5) ){
  axis( 2, cex.axis=1.5)
  mtext( line=4, side=2, "Ystuff")}
}

# same trick of adding a legend strip. 
par( oma=c(0,0,0,0))
image.plot( zlim=zr, legend.only=TRUE, horizontal=TRUE, legend.mar=5)

# reset panel 
set.panel()


####
# show colors
## the factory colors:

clab<- colors()
n<- length( clab)
N<- ceiling( sqrt(n) )
M<- N
temp<- rep( NA,M*N)
temp[1:n] <- 1:n
z<- matrix(temp, M,N)

image(seq(.5,M+.5,,M+1), seq(.5,N+.5,,N+1)
       , z,  col=clab, axes=FALSE, xlab="", ylab="")


# see the  function  fields.color.picker() to locate colors



# dumping out colors by name for a latex document
# this creates text strings that are the LaTeX color definitions
# using the definecolor function. 

# grab all of the R default colors 
clab<- colors()

for( nn in clab){
  temp<- signif(col2rgb(nn)/256, 3)
   cat(
    "\definecolor{", 
                nn, "}",
    "{rgb}{", temp[1], 
          ",", temp[2], 
          ",", temp[3], 
           "}", fill=TRUE , sep="")
 }

# this loop prints out definitions such as 
# \definecolor{yellowgreen}{rgb}{0.602,0.801,0.195}
# having loaded the color package in LaTeX 
# defining this color 
# use the construction  {\color{yellowgreen}  THIS IS A COLOR}
# to use this color in a talk or document. 

# this loop prints out all the colors in LaTeX language
# as their names and can be converted to a pdf for handy reference.

sink( "showcolors.tex")

clab<- colors()
for( nn in clab){
  temp<- signif(col2rgb(nn)/256, 3)
   cat(
    "\definecolor{", 
                nn, "}",
    "{rgb}{", temp[1], 
          ",", temp[2], 
          ",", temp[3], 
           "}", fill=TRUE , sep="")
   cat( paste("{ \color{",nn,"} ", nn," $\bullet$ \\ }", sep=""),
                    fill=TRUE)
}
sink()

Description

The characteristics of an ionizing flame are varied with the intent of maximizing the intensity of emitted light for lithuim in solution. Areas outside of the measurements are where the mixture may explode! Note that the optimum is close to the boundary. Source of data is from a master's level lab experiment in analytical chemistry from Chuck Boss's course at NCSU. <s-section name= "DATA DESCRIPTION"> This is list with the following components

Arguments

Description

A moderate size (about 8400 locations) spatial dataset that is well-known in the applied mathematics approximation literature for testing interpolation methods.

Usage

data(glacier)

Format

The format of glacier is a list with two components:

loc:

8338x2 matrix of the locations (meters??).

y:

A vector of elevations (meters ??).

Details

This data set appears in papers that develop interpolation methods for scattered data and serves as an interesting bridge to the examples in applied math that develop radial basis function surface fitting. The data was originally used by R. Franke.

Unfortunately at this time we can not find any background on where these data were collected or indeed even the location of this glacier. However, it is an interesting data set in that it appears that the elevations are reported along lnes of equal elevation, i.e. contours, perhaps from a digitization of a topographic map or survey. It is important to estimate the surface in a way that the artifacts from discretization are not present. In the example below the compactly supported kernel interpolation still has some artifacts.

The glacier data set is available at this website https://oleg-davydov.de/scat_data.html

The examples below are useful for comparing different approximations to a Gaussian spatial process estimate for the elevation surface. Of course in using a stationary covariance ( e.g. the Matern or Wendland) these are also radial basis smoothing or interpolation of the data.

Examples

data( glacier )
# EDA for raw obs:

bubblePlot( glacier$loc, glacier$y, highlight=FALSE, size=.5)

# identifying contour levels. Note this is reported at regular levels
# (Every 25m ???)

table( glacier$y)



# find sigma and rho by maximum likelihood 
# for a fixed range
#  the default is the Wendland covariance with k=2
# See help(Wendland)

# this takes about 5 minutes
# macbook pro Quad-Core Intel Core i5 8 GB

#options(spam.nearestdistnnz=c(5e7,1e3))
#system.time( 
# obj0<- fastTps(glacier$loc, glacier$y, 
#                       theta=2,
#                      profileLambda=TRUE) 
#)
 

# set.panel(2,2)
# plot( obj0)
# set.panel()

# just evaluate at MLE
# reset default matrix size that the spam pacakge will use.

## Not run: 

options(spam.nearestdistnnz=c(5e7,1e3))
system.time( obj1<- 
               fastTps(glacier$loc, glacier$y, 
                       theta=2,
                       lambda= 7.58e-5
                        ) 
)

system.time(
look1<- predictSurface( obj1, nx=150, ny=150)
)

imagePlot( look1)


system.time(
out<- simLocal.spatialProcess(obj1, M=3, nx=150, ny=150)
)
set.panel( 2,2)
imagePlot( look1)
zlim<- range( out$z, na.rm=TRUE)
for( k in 1:3){
imagePlot(out$x, out$y, out$z[,,k], zlim=zlim)
}

# near interpolation surface using Matern smoothness .5 
 system.time( 
 obj2<- spatialProcess(glacier$loc, glacier$y,
                          aRange = 1.5, 
                          lambda = 1e-5,
                          smoothness = .5)
 )
 
system.time(
out<- simLocal.spatialProcess(obj2, M=3, nx=150, ny=150,
fast=TRUE)
)

set.panel( 2,2)
imagePlot( look1)
zlim<- range( out$z, na.rm=TRUE)
for( k in 1:3){
imagePlot(out$x, out$y, out$z[,,k], zlim=zlim)
}

system.time(
look2<- predictSurface.mKrig( obj2, nx=150, ny=150,
                fast=TRUE, NNSize=5)
)

system.time(
look2B<- predictSurface( obj2, nx=150, ny=150,
                fast=FALSE)
)

err<- c((look2$z - look2B$z)/look2B$z)
stats( log10( abs(err) ) )

# some error plots ( percent relative error)
imagePlot(look2$x, look2$y, 100*(look2$z - look2B$z)/look2B$z  )

imagePlot(look2$x, look2$y, 100*(look1$z - look2B$z)/look2B$z  )


## End(Not run)

Description

The object grid.list refers to a list that contains information for evaluating a function on a 2-dimensional grid of points. If a function has more than two independent variables then one also needs to specify the constant levels for the variables that are not being varied. This format is used in several places in fields for functions that evaluate function estimates and plot surfaces. These functions provide some default conversions among information and the gird.list. The function discretize.image is a useful tool for "registering" irregular 2-d points to a grid.

Usage

makeMultiIndex(M)
parse.grid.list( grid.list, order.variables="xy")
fields.x.to.grid(x,nx=80, ny=80, xy=c(1,2))
fields.convert.grid( midpoint.grid )
discretize.image(x, m = 64, n = 64, grid = NULL, 
                 expand = c(1 + 1e-08, 1 + 1e-08),
                 boundary.grid = FALSE, na.rm = TRUE)
make.surface.grid( grid.list)
unrollZGrid( grid.list, ZGrid)

Arguments

Details

makeMultiIndex creates an expanded set of indices to referencce a regular grid. M are L integers with product prodM Will create a prodM by L matrix that is all combinations of (1:M[i]) for i =1,2, ...L This is organized in the standard array ordering where the first column varies the fastest for M = c( 3,2,4) the result will be a 24X3 matrix with the entries:

                          1,1,1 
                          2,1,1
                          3,1,1
                          1,2,1
                          2,2,1
                          3,2,1
                      etc ...
		          ...
             and  ending with
                          2,2,4
                          3,2,4
 

All about grid lists:

The form of a grid.list is

list( var.name1= what1 , var.name2=what2 , ... var.nameN=what3)

Here var.names are the names of the independent variables. The what options describe what should be done with this variable when generating the grid. These should either an increasing sequence of points or a single vaules. Obviously there should be only be two variables with sequences to define a grid for a surface.

Most of time the gridding sequences are equally spaced and are easily generated using the seq function. Also throughout fields the grid points are typically the midpoints of the grid rather the grid box boundaries. However, these functions can handle unequally spaced grids and the logical boundary.grid can indicate a grid being the box boundaries.

The variables in the list components are assumed to be in the same order as they appear in the data matrix.

A useful function that expands the grid from the grid.list description into a full set of locations is make.surface.grid and is just a wrapper around the R base function expand.grid. A typical operation is to go from a grid.list to the set of grid locations. Evaluate a fucntion at these lcoations and then reformat this as an image for plotting. Here is how to do this cleanly:

grid.list<- list( x= 1:10, y=1:15)
xg<- make.surface.grid(grid.list)
# look at a surface dependin on xg locations
z<- xg[,1] + 2*xg[,2]
out<- list( x=grid.list$x, y= grid.list$y, z=matrix( z, nrow=10, ncol=15))
# now for example
image.plot( out)

The key here is that xg and matrix both organize the grid in the same order.

Some fields internal functions that support interpreting grid list format are:

fields.x.to.grid: Takes an "x" matrix of locations or independent variables and creates a reasonable grid list. This is used to evaluate predicted surfaces when a grid list is not explicited given to predictSurface. The variables (i.e. columns of x) that are not part of the grid are set to the median values. The x grid values are nx equally spaced points in the range x[, xy[1]]. The y grid values are ny equally spaced points in the range x[, xy[2]].

parse.grid.list: Takes a grid list and returns the information in a more expanded list form that is easy to use. This is used, for example, by predictSurface to figure out what to do!

fields.convert.grid: Takes a vector of n values assumed to be midpoints of a grid and returns the n+1 boundaries. See how this is used in discretize.image with the cut function. This function will handle unequally spaced grid values.

discretize.image: Takes a vector of locations and a 2-d grid and figures out to which boxes they belong. The output matrix ind has the grid locations. If boundary.grid is FALSE then the grid list (grid) is assumed to be grid midpoints. The grid boundaries are taken to be the point half way between these midpoints. The first and last boundaries points are determined by extrapolating so that the first and last box has the midpoint in its center. (See the code in fields.convert.grid for details.) If grid is NULL then midpoints are found from m and n and the range of the x matrix.

unrollZGrid Checks that the ZGrid object is compatible with th e grid.list and concatenates the grid arrays into vectors. This version of the covariates are used the usual predict function.

See Also

as.surface, predictSurface, plot.surface, surface, expand.grid, as.image

Examples

#Given below are some examples of grid.list objects and the results  
#when they are used with make.surface.grid. Note that  
#make.surface.grid returns a matrix that retains the grid.list  
#information as an attribute. 

grid.l<- list( 1:3, 2:5) 
make.surface.grid(grid.l)

  
grid.l <- list( 1:3, 10, 1:3) 
make.surface.grid(grid.l) 

#The next  example shows how the grid.list can be used to  
#control surface plotting and evaluation of an estimated function. 
# first create a test function  

set.seed( 124)

X<- 2*cbind( runif(30), runif(30), runif(30)) -1
  
dimnames( X)<- list(NULL, c("X1","X2","X3")) 
y<- X[,1]**2 + X[,2]**2 + exp(X[,3])   

# fit an  interpolating thin plate spline  
out<- Tps( X,y) 

grid.l<- list( X1= seq( 0,1,,20), X2=.5, X3=seq(0,1,,25)) 
surface( out, grid.list=grid.l) 
#  surface plot based on a 20X25 grid in X1 an X3  
#                       over the square [0,2] and [0,2]   
#                       holding X2 equal to 1.0. 
#


# test of discretize to make sure points on boundaries are counted right
set.seed(123)
x<- matrix( runif(200), 100,2)
look<- discretize.image( x, m=2,n=2)
xc<- seq(min(x[,1]), max(x[,1]),,5)
xc<- xc[2:4]
yc<- seq(min(x[,2]), max(x[,2]),,5)
yc<- yc[2:4]
grid <-  list( x= xc, y= yc)
look2<- discretize.image( x, m=2,n=2)

table( look$index )
table( look2$index )

# indicator image of discretized locations
look<- discretize.image( RMprecip$x, m=15, n=15)
image.plot( look$grid$x, look$grid$y,look$hist )  
# actual locations
points( RMprecip$x,col="magenta", pch=".")

Description

Given two sets of locations defined on a 2-d grid efficiently multiplies a cross covariance with a vector. The intermediate compuations (the setup) can also be used for fast simulation of the processes on a grid using the circulant embedding technique.

Usage

stationary.image.cov(ind1, ind2, Y, cov.obj = NULL, setup = FALSE, 
grid, M=NULL,N=NULL,cov.function="stationary.cov", delta = NULL, cov.args = NULL, ...) 

Exp.image.cov(ind1, ind2, Y, cov.obj = NULL, setup = FALSE, grid, ...)

Rad.image.cov(ind1, ind2, Y, cov.obj = NULL, setup = FALSE, grid, ...)

matern.image.cov(ind1, ind2, Y, cov.obj = NULL, setup = FALSE, grid,
M=NULL,N=NULL,aRange= 1.0, smoothness=.5, theta=NULL)

wendland.image.cov(ind1, ind2, Y, cov.obj = NULL, 
    setup = FALSE, grid, M = NULL, N = NULL, cov.args=NULL, ...)

Arguments

Details

This function was provided to do fast computations for large numbers of spatial locations and supports the conjugate gradient solution in krig.image. In doing so the observations can be irregular spaced but their coordinates must be 2-dimensional and be restricted to grid points. (The function as.image will take irregular, continuous coordinates and overlay a grid on them.)

Returned value: If ind1 and ind2 are matrices where nrow(ind1)=m and nrow(ind2)=n then the cross covariance matrix, Sigma, is an mXn matrix (i,j) element is the covariance between the grid locations indexed at ind1[i,] and ind2[j,]. The returned result is Sigma multiplied by Y. Note that one can always recover the coordinates themselves by evaluating the grid list at the indices. E.g. If x and y are the grids for the X and Y dimensions, cbind( x[ind1[,1]], y[ind1[,2])) will give the coordinates associated with ind1. Clearly it is better just to work with ind1!

Functional Form: Following the same form as Exp.cov stationary.cov for irregular locations, the covariance is defined as phi( D.ij) where D.ij is the Euclidean distance between x1[i,] and x2[j,] but having first been scaled by aRange. Specifically,

D.ij = sqrt( sum.k (( x1[i,k] - x2[j,k]) /aRange[k])**2 ).

See Matern for the version of phi for the Matern family.

Note that if aRange is a scalar then this defines an isotropic covariance function.

Implementation: This function does the multiplication on the full grid efficiently by a 2-d FFT. The irregular pattern in Y is handled by padding with zeroes and once that multiplication is done only the appropriate subset is returned.

As an example assume that the grid is 100X100 let big.Sigma denote the big covariance matrix among all grid points ( If the parent grid is 100x100 then big.Sigma is 10K by 10K !) Here are the computing steps:

temp<- matrix( 0, 100,100)

temp[ ind2] <- Y

temp2<- big.Sigma%*% temp

temp2[ind1]

Notice how much we pad with zeroes or at the end throw away! Here the matrix multiplication is effected through convolution/FFT tricks to avoid creating and multiplying big.Sigma explicitly. It is often faster to multiply the regular grid and throw away the parts we do not need then to deal directly with the irregular set of locations.

Note: In this entire discussion Y is treated as vector. However if one has complete data then Y can also be interpreted as a image matrix conformed to correspond to spatial locations. See the last example for this distinction.

Value

A vector that is the multiplication of the cross covariance matrix with the vector Y.

See Also

smooth.2d, as.image, krig.image, stationary.cov

Examples

# multiply 2-d isotropic exponential with aRange=4 by a random vector 

junk<- matrix(rnorm(100*100), 100,100)

cov.obj<- stationary.image.cov( setup=TRUE, 
             grid=list(x=1:100,y=1:100),aRange=8) 
result<-  stationary.image.cov(Y=junk,cov.obj=cov.obj)

image( matrix( result, 100,100)) # NOTE that is also a smoother!

# to do it again, no setup is needed 
#  e.g. 
#  junk2<- matrix(rnorm(100**2, 100,100))
#  result2<-  stationary.image.cov(Y=junk2, cov.obj=cov.obj)

# generate a grid and set of indices based on discretizing the locations
# in the precip dataset

 out<-as.image( RMprecip$y, x= RMprecip$x)
 ind1<- out$ind
 grid<- list( x= out$x, y=out$y)

#
# discretized x locations  to use for comparison
  xd<- cbind( out$x[ out$ind[,1]], out$y[ out$ind[,2]] )

# setup to create cov.obj for exponential covariance with range= 1.25

 cov.obj<- stationary.image.cov( setup=TRUE, grid=grid, aRange=1.25) 

# multiply covariance matrix by an arbitrary vector
 junk<-  rnorm(nrow( ind1))
 result<- stationary.image.cov( ind1, ind1, Y= junk,cov.obj=cov.obj)

# The brute force way would be  
#   result<- stationary.cov( xd, xd, aRange=1.25, C=junk)
# or 
#   result<- stationary.cov( xd, xd, aRange=1.25) %*% junk
# both of these take much longer 


# evaluate the covariance between all grid points and the center grid point
 Y<- matrix(0,cov.obj$m, cov.obj$n)
 Y[32,32]<- 1
 result<- stationary.image.cov( Y= Y,cov.obj=cov.obj)
# covariance surface with respect to the grid point at (32,32)
# 
# reshape "vector" as an image
 temp<-  matrix( result, cov.obj$m,cov.obj$n)
 image.plot(cov.obj$grid$x,cov.obj$grid$y, temp)
# or persp( cov.obj$grid$x,cov.obj$grid$y, temp) 

# check out the Matern
 grid<- list( x= seq(-105,-99,,64), y= seq( 40,45,,64)) 
 cov.obj<- matern.image.cov( 
             setup=TRUE, grid=grid, aRange=.55, smoothness=1.0)
 Y<- matrix(0,64,64)
 Y[16,16]<- 1

 result<- matern.image.cov( Y= Y,cov.obj=cov.obj)
  temp<-  matrix( result, cov.obj$m,cov.obj$n)
 image.plot( cov.obj$grid$x,cov.obj$grid$y, temp)

# Note we have centered at the location (grid$x[16],grid$y[16]) for this case
#  using sim.rf to simulate an Matern field
  look<- sim.rf( cov.obj)
  image.plot( grid$x, grid$y, look)

Description

This function combines the R image function with some automatic placement of a legend. This is done by splitting the plotting region into two parts. Putting the image in one and the legend in the other. After the legend is added the plot region is reset to the image plot. This function also allows for plotting quadrilateral cells in the image format that often arise from regular grids transformed with a map projection or a scaling and rotation of coordinates. See the example where this function can create a similar graphic to the ggplot package. image.plot functionality has been frozen, see the more recent function imagePlot which is backwardly compatible with this function.

Usage

## S3 method for class 'plot'
image(...,
                 add = FALSE, 
                 breaks= NULL, nlevel = 64, col = NULL,  
    horizontal = FALSE, legend.shrink = 0.9, legend.width = 1.2, 
    legend.mar = ifelse(horizontal, 3.1, 5.1), legend.lab = NULL,
    legend.line= 2,                    
    graphics.reset = FALSE, bigplot = NULL, smallplot = NULL, 
    legend.only = FALSE,  lab.breaks = NULL, 
    axis.args = NULL, legend.args = NULL, legend.cex=1.0,
    midpoint = FALSE, border = NA, 
    lwd = 1,verbose = FALSE )

Arguments

Details

This is a function using the basic R graphics. The coding was done to make it easier for users to see how this function works and to modify.

How this function works: The strategy for image.plot is simple, divide the plotting region into two smaller regions bigplot and smallplot. The image goes in one and the legend in the other. This way there is always room for the legend. Some adjustments are made to this rule by not shrinking the bigplot if there is already room for the legend strip and also sticking the legend strip close to the image plot. One can specify the plot regions explicitly by bigplot and smallplot if the default choices do not work.(Note that these in figure coordinates. ) There may be problems with small plotting regions in fitting both of these elements into the plot region and one may have to change the default character sizes or margins to make things fit. Sometimes this function will not reset the type of margins correctly and the "null" call par(mar = par("mar")) may help to fix this issue.

The text is too small! This always seems to happen as one is rushing to finish a talk and the figures have tiny default axis labels. Try just calling the function fields.style before plotting. List out this function to see what is changed, however, all text is increased by 20% in size.

Why “image.plot" and not “image"? The R Base function image is very useful but it is awkward to place a legend quickly. However, that said if you are drawing several image plots and want a common legend use the image function and just just use image.plot to add the legend. See the example in the help file. Note that you can use image to draw a bunch of images and then follow with image.plot and legend.only=TRUE to add a common legend. (See examples below.)

Almost cloropleths too: It should be noted that this image function is slightly different than a cloropleth map because the legend is assuming that a continous scale has been discretized into a series of colors. To make the image.plot function as a cloropleth graphic one would of course use the breaks option and for clarity perhaps code the different regions as different integer values. In addition, for publication quality one would want to use the legend.args to add more descriptive labels at the midpoints in the color strip.

Relationship of x, y and z: If the z component is a matrix then the user should be aware that this function locates the matrix element z[i,j] at the grid locations (x[i], y[j]) this is very different than simply listing out the matrix in the usual row column tabular form. See the example below for details on the difference in formatting. What does one do if you do not really have the "z" values on a regular grid? See the functions quilt.plot.Rd and as.image to discretise irregular observations to a grid. If the values makes sense as points on a smooth surface see Tps and fastTps for surface interpolation.

Adding separate color to indicate the grid box boundaries. Sometimes you want to show to the grid box borders to emphasize this is not a smooth surface. To our knowledge there is no easy way to do this as an option in image. But if your image is formatted in the "poly image" style where x and y are also matrices you can use the polyimage (see the border argument above) option to draw in boundaries.

Grids with unequally spacing – quadrialteral pixels: If x and y are matrices that are a smooth transformation of a regular grid then z[i,j] can be interperted as representing the average value in a quadrilateral that is centered at x[i,j] and y[i,j] (midpoint TRUE). The details of how this cell is found are buried in poly.image but it it essentially found using midpoints between the centers. If midpoint is FALSE then x and y are interpreted as the corners of the quadrilateral cells. But what about z? The four values of z are now averaged to represent a value at the midpoint of the cell and this is what is used for plotting. Quadrilateral grids were added to help with plotting the gridded output of geophysical models where the regular grid is defined according to one map projection but the image plotting is required in another projection. Typically the regular grid becomes distorted in a smooth way when this happens. See the regional climate example for a illustration of this application. One can add border colors in this case easily because these choices are just passed onto the polygon function.

Adding the pixel grid for rectangular images: For adding the grid of pixel borders to a rectangular image try this example after calling image.plot.

 
 dx <- x[2] - x[1]  
 dy <- y[2] - y[1]  
 xtemp<- seq(  min( x)- dx/2, max(x)+ dx/2,
         length.out = length(x) +1) 
 ytemp<- seq(  min( y)- dy/2, max(y)+ dy/2,
         length.out = length(y) +1)
 xline( xtemp, col="grey", lwd=2)
 yline( ytemp, col="grey", lwd=2) 

Here x and y here are the x and y grid values from the image list.

Fine tuning color scales: This function gives some flexibility in tuning the color scale to fit the rendering of z values. This can either be specially designed color scale with specific colors ( see help on designer.colors), positioning the colors at specific points on the [0,1] scale, or mapping distinct colors to intervals of z. The examples below show how to do each of these. In addition, by supplying lab.break strings or axis parameters one can annotate the legend axis in an informative matter.

Adding just the legend strip: Note that to add just the legend strip all the numerical information one needs is the zlim argument and the color table! See examples for tricks in positioning.

About color tables: We like tim.colors as a default color scale and so if this what you use this can be omitted. Unfortunately this is not the default for the image function. Another important color scale is viridis() from the viridis package. It seems that by and large everyone seems to react postively to viridis – guess that is the point!

The topographic color scale (topo.colors) is also a close second showing our geophysical bias. Users may find larry.colors useful for coding distinct regions in the style of a cloropleith map. See also terrain.colors for a subset of the topo ones and designer.colors to "roll your own" color table. One nice option in this last function is to fix color transitions at particular quantiles of the data rather than at equally spaced intervals. For color choices see how the nlevels argument figures into the legend and main plot number of colors. Also see the colors function for a listing of all the colors that come with the R base environment.

The details of placing the legend and dividing up the plotting real estate: It is surprising how hard it is to automatically add the legend! All "plotting coordinates" mentioned here are in device coordinates. The plot region is assumed to be [0,1]X[0,1] and plotting regions are defined as rectangles within this square. We found these easier to work with than user coordinates.

legend.width and legend.mar are in units of character spaces. These units are helpful in thinking about axis labels that will be put into these areas. To add more or less space between the legend and the image plot alter the mar parameters. The default mar settings (5.1,5.1,5.1,2.1) leaves 2.1 spaces for vertical legends and 5.1 spaces for horizontal legends.

There are always problems with default solutions to placing information on graphs but the choices made here may be useful for most cases. The most annoying thing is that after using image.plot and adding information the next plot that is made may have the slightly smaller plotting region set by the image plotting. The user should set reset.graphics=TRUE to avoid the plotting size from changing. The disadvantage, however, of resetting the graphics is that one can no longer add additional graphics elements to the image plot. Note that filled.contour always resets the graphics but provides another mechanism to pass through plotting commands. Apparently filled.contour, while very pretty, does not work for multiple plots.

About setup and add legend functions These came about to create a scatterplot in Base R Graphics where the points are colored with a color scale and the scale can be plotted as part of the figure See bubblePlot for a version of this kind of figure. The function setupLegend should be used first to create enough space to add a color scale later. After plotting then addLegend will add the color scale. Note that if the color scale has been created by the color.scale function the last call to this function will use the color scale and limits created in color.scale.

In summary here is an example of using these functions with the colors in mind:

  info<- setupLegend()
  colTab<- rainbow(10)
  plot( 1:10, 201:210, col=colTab, pch=16)
  addLegend(info, col=colTab, zlim = c(1,10))

Here is one where four colors are mapped to specific values (ala image).

  info<-setupLegend()
  colTab= color.scale(201:210, rainbow(4))
  plot( 1:10, 201:210, col=colTab, pch=16 )
  addLegend(info, col=colTab, zlim = c(201,210) )

More complete graphics languages, such as that in ggplot, do not need such functions because the entire graphics segment is parsed to create the complete figure. In this way room for a color scale can be created automatically. The functions proposed here are a simple work around to create these figures using base R graphics.

Other packages levelplot that is part of the lattice package has a very similar function to image.plot and a formula syntax in the call. The geom_raster for setting up a graphics object within ggplot is another alternative forr image plots with legends. See the last example to compare the steps in creating an image plot using image.plot that is close to the ggplot version. Mostly this involves resetting base graphics parameters using the par function.

Multiple images: By keeping the zlim argument the same across images one can generate the same color scale. (See the image help file.) One useful technique for a panel of images is to just draw the images with good old image and then use image.plot to add a legend to the last plot. (See example below for messing with the outer margins to make this work.) Usually a square plot (pty="s") done in a rectangular plot region will have room for the legend stuck to the right side without any other adjustments. See the examples below for more complicated arrangements of multiple image plots and a summary legend. The reader is also referred to the package autoimage as a set of functions in base to help with drawing multiple images and also more support for geographic coordinates.

Side Effects

After exiting, the plotting region may be changed to make it possible to add more features to the plot. To be explicit, par()\$plt may be changed to reflect a smaller plotting region that has accommodated room for the legend subplot.

If xlim and ylim are specified the pixels may overplot the axis lines. Just use the box function to redraw them.

See Also

imagePlot, image,poly.image, filled.contour, quilt.plot, bubblePlot, plot.surface, add.image, colorBar, tim.colors, designer.colors

Examples

x<- 1:10
  y<- 1:15
  z<- outer( x,y,"+") 
  image.plot(x,y,z) 

# or 
  obj<- list( x=x,y=y,z=z)
  image.plot(obj, legend.lab="Sverdrups")
  
################################################################ 
# the next sequence of examples explain how to quickly 
# adpat this basic plot to include morre features
# In another direction see the very last example where 
# we use many of the setting in base R graphic to mimic a 
# (beautiful) ggplot version. 
###############################################################
#
# add some points on diagonal using standard plot function
#(with some clipping beyond 10 anticipated)

  points( 5:12, 5:12, pch="X", cex=3)
  
# in general image.plot will reset the plot window so you
# can add any feature that normally works in base R
# e.g. lines, text, contour, boxplots, ....
#
# adding breaks and distinct colors for intervals of z
# with and without lab.breaks

  brk<- quantile( c(z))
  image.plot(x,y,z, breaks=brk, col=rainbow(4))
  
# annotate legend strip with the  break point values and add a label

  image.plot(x,y,z, breaks=brk, col=rainbow(4),
                       lab.breaks=names(brk))
#
# compare to 

  zp <-quantile(c(z), c( .05, .1,.5, .9,.95))
  image.plot(x,y,z, 
     axis.args=list( at=zp, labels=names(zp) ) )
     
# a log scaling for the colors

  ticks<- c( 1, 2,4,8,16,32)
  image.plot(x,y,log(z), axis.args=list( at=log(ticks), labels=ticks))

# see help file for designer.colors to generate a color scale that adapts to 
# quantiles of z. 
# Add some color scales together here is an example of  5 blues to white to 5 reds
# with white being a specific size.
 colorTable<- designer.colors(11, c( "blue","white", "red") )
# breaks with a gap of 10 to 17 assigned the white color
 brks<- c(seq( 1, 10,,6), seq( 17, 25,,6)) 
 image.plot( x,y,z,breaks=brks, col=colorTable)
#
#fat (5 characters wide) and short (50% of figure)  color bar on the bottom
   image.plot( x,y,z,legend.width=5, legend.shrink=.5, horizontal=TRUE) 

# adding a label with all kinds of additional arguments.
# use side=4 for vertical legend and side= 1 for horizontal legend
# to be parallel to axes. See help(mtext).

image.plot(x,y,z, 
       legend.args=list( text="unknown units",
     col="magenta", cex=1.5, side=4, line=2))
     
# and finally add some grid lines
 dx <- x[2] - x[1]  
 dy <- y[2] - y[1]  
 xtemp<- seq(  min( x)- dx/2, max(x)+ dx/2,
         length.out = length(x) +1) 
 ytemp<- seq(  min( y)- dy/2, max(y)+ dy/2,
         length.out = length(y) +1)
 xline( xtemp, col="grey", lwd=2)
 yline( ytemp, col="grey", lwd=2)

###############################################################
#### example using an irregular quadrilateral grid
###############################################################
data( RCMexample)

image.plot( RCMexample$x, RCMexample$y, RCMexample$z[,,1])
ind<- 50:75 # make a smaller image to show bordering lines
image.plot( RCMexample$x[ind,ind], RCMexample$y[ind,ind], RCMexample$z[ind,ind,1],
                                      border="grey50", lwd=2)

###############################################################
#### multiple images with a common legend
###############################################################
set.panel()

# Here is quick but quirky way to add a common legend to several plots.
# The idea is leave some room in the margin and then at the end 
# overplot the legend in this margin

par(oma=c( 0,0,0,4)) # margin of 4 spaces width at right hand side
set.panel( 2,2) # 2X2 matrix of plots

# now draw all your plots using usual image command
for (  k in 1:4){
  data<- matrix( rnorm(150), 10,15)
  image( data, zlim=c(-4,4), col=tim.colors())
# and just for fun add a contour plot  
  contour( data, add=TRUE)
}

par(oma=c( 0,0,0,1))# reset margin to be much smaller.
image.plot( legend.only=TRUE, zlim=c(-4,4)) 

# image.plot tricked into  plotting in margin of old setting 

set.panel() # reset plotting device

#
# Here is a more learned strategy to add a common legend to a panel of
# plots  consult the split.screen help file for more explanations.
# For this example we draw two
# images top and bottom and add a single legend color bar on the right side 

# first divide screen into the figure region (left) and legend region (right)
   split.screen( rbind(c(0, .8,0,1), c(.8,1,0,1)))

# now subdivide up the figure region into two parts
   split.screen(c(2,1), screen=1)-> ind
   zr<- range( 2,35)
# first image
   screen( ind[1])
   image( x,y,z, col=tim.colors(), zlim=zr)

# second image
   screen( ind[2])
   image( x,y,z+10, col=tim.colors(), zlim =zr)

# move to skinny region on right and draw the legend strip 
   screen( 2)
   image.plot( zlim=zr,legend.only=TRUE, smallplot=c(.1,.2, .3,.7),
   col=tim.colors())

   close.screen( all=TRUE)


# you can always add a legend arbitrarily to any plot;
# note that here the plot is too big for the vertical strip but the
# horizontal fits nicely.
plot( 1:10, 1:10)
image.plot( zlim=c(0,25), legend.only=TRUE)
image.plot( zlim=c(0,25), legend.only=TRUE, horizontal =TRUE)

# combining the  usual image function and adding a legend
# first change margin for some more room
## Not run: 
par( mar=c(10,5,5,5))
image( x,y,z, col=topo.colors(64))
image.plot( zlim=c(0,25), nlevel=64,legend.only=TRUE, horizontal=TRUE,
col=topo.colors(64))

## End(Not run)
#

# adding a legend by  automatically making room. 
# and coloring points
  info<- setupLegend()
  colTab<- rainbow(10)
  plot( 201:210, 201:210, col=colTab, pch=16)
  addLegend(info, col=colTab, zlim = c(201,210))
#

#######################################################
##### Comparison to ggplot
#######################################################
# the following example was created as way avoid doing more important
# things 
# Note how close base graphics can get to reproducing the ggplot style.

## Not run: 
library( viridis)
library(ggplot2)

x<- 1:20
y<-  1:24
z<- outer( x, y, "+")


# ggplot version 
  mesh<- expand.grid( x= x, y=y)
  mesh$z <- c(z)
  ggplot( data=mesh, aes( x=x, y=y, fill=z)) + 
    geom_raster(interpolate= FALSE)  + 
    scale_fill_continuous(type = "viridis")  +
    theme_bw()


# inflate range to give a margin around image
  xr<- range(x) +  c(-.08, .08)* diff( range(x))
  yr<- range(y) +  c(-.08, .08)* diff( range(y))
  
# changing these graphics parameters tends to push 
# text closer to the axes. 
  par( mgp=c(1.5,.5,0),mar=c(2.5,2.5,.5,1), cex=.8)
  
  image.plot(x,y,z, 
             col = viridis(128), 
   legend.shrink = .27,
            xlim = xr, 
            ylim = yr,
    legend.width = 1.5,
      legend.mar = 3,
     legend.args = list( text = "z",
                          cex = .8,
                         side = 3,
                         line = .5)
   )


## End(Not run)