mapImage: Plot an image on a existing map

Description

Plot an image on an existing map

Usage

mapImage(longitude, latitude, z, zlim, zclip=FALSE,
         breaks, col, colormap, border=NA,
         lwd=par("lwd"), lty=par("lty"),
         filledContour=FALSE, missingColor=NA, debug=getOption("oceDebug"))

Arguments

longitude

longitudes of grid lines

latitude

latitudes of grid lines

matrix for image

zlim

limit for z (colour)

zclip

logical, indicating whether to clip the colours by using the missingColor for out-of-range values. If if zlim is given then its min and max set the range. If zlim is not given but brea

breaks

the z values for breaks in the colour scheme. If this is of length 1, the value indicates the desired number of breaks, which is supplied to pretty, in determining clean break points.

col

either a vector of colours corresponding to the breaks, of length 1 plus the number of breaks, or a function specifying colours, e.g. oce.colorsJet for a rainbow.

colormap

optional colormap, as created by colormap. If a colormap is provided, then its properties takes precedence over breaks, col, missingColor

border

colour used for borders of patches (passed to polygon); the default NA means no border.

lwd

line width, used if borders are drawn

lty

line type, used if borders are drawn

filledContour

either a boolean value indicating whether to use filled contours to plot the image, or a numerical value indicating the resampling rate to be used in interpolating from lon-lat coordinates to x-y coordinates. See Details

missingColor

a color to be used to indicate missing data, or NA to skip the drawing of such regions (which will retain whatever material has already been drawn at the regions).

debug

a flag that turns on debugging. Set to 1 to get a moderate amount of debugging information, or to 2 to get more. (Temporary note: setting this to 99 forces the use of a drawing method that draws polygons one by one, which is possi

Details

Adds an image to an existing map, by analogy to image.

The data are on a regular grid in lon-lat space, but not in the projected x-y space. This means that image cannot be used. Instead, there are two approaches, depending on the value of filledContour.

If filledContour is FALSE, the image ``pixels'' are with polygon. This can be slow for fine grids. Also, the results can seem ``washed out'' on some plotting devices, owing to an issue with anti-aliasing in the drawing of the pixel polygons. If a png device is to be used, it is advised to supply arguments type="cairo" and antialias="none"; see [1].

If filledContour is TRUE or a numerical value, then the the ``pixels'' are remapped into a regular grid and then displayed with .filled.contour. The remapping starts by converting the regular lon-lat grid to an irregular x-y grid using lonlat2map. This irregular grid is then interpolated onto a regular x-y grid with interp from the akima package. If filledContour is TRUE, the dimensions of the regular x-y grid is the same as that of the original lon-lat grid; otherwise, the number of rows and columns are multiplied by the numerical value of filledContour, e.g. the value 2 means to make the grid twice as fine. Filling contours can produce aesthetically-pleasing results, but the method involves interpolation, so the data are not represented exactly and analysts are advised to compare the results from the two methods (and perhaps various grid refinement values) to guard against misinterpretation.

References

1. http://codedocean.wordpress.com/2014/02/03/anti-aliasing-and-image-plots/

Examples

Run this code

library(oce)
data(coastlineWorld)
data(topoWorld)

par(mfrow=c(2,1), mar=c(2, 2, 1, 1))
lonlim <- c(-70,-50)
latlim <- c(40,50)
topo <- decimate(topoWorld, by=2) # coarse to illustrate filled contours
topo <- subset(topo, latlim[1] < latitude & latitude < latlim[2])
topo <- subset(topo, (360+lonlim[1]) < longitude & longitude < (360+lonlim[2]))
mapPlot(coastlineWorld, type='l',
        longitudelim=lonlim, latitudelim=latlim,
        proj="polyconic", orientation=c(90,-60,0), grid=TRUE)
breaks <- seq(-5000, 1000, 500)
mapImage(topo, col=oce.colorsGebco, breaks=breaks)
mapLines(coastlineWorld)
box()
mapPlot(coastlineWorld, type='l',
        longitudelim=lonlim, latitudelim=latlim,
        proj="polyconic", orientation=c(90,-60,0), grid=TRUE)
mapImage(topo, filledContour=TRUE, col=oce.colorsGebco, breaks=breaks)
box()
mapLines(coastlineWorld)

## Northern polar region, with colour-coded bathymetry
drawPalette(c(-5000,0), zlim=c(-5000, 0), col=oce.colorsJet)
mapPlot(coastlineWorld, type='l',
        longitudelim=c(-180,180), latitudelim=c(60,120),
        proj="stereographic", grid=TRUE)
mapImage(topoWorld, zlim=c(-5000, 0), col=oce.colorsJet)
mapLines(coastlineWorld[['longitude']], coastlineWorld[['latitude']])

# Levitus SST (requires dataset from http://www.esrl.noaa.gov)
par(mfrow=c(1,1))
library(ncdf4)
con <- nc_open("/data/oar/levitus/temperature_annual_5deg.nc")
lon <- ncvar_get(con, "lon")
lat <- ncvar_get(con, "lat")
SST <- ncvar_get(con, "t_mn")[,,1]
par(mar=rep(1, 4))
Tlim <- c(-2, 30)
drawPalette(Tlim, col=oce.colorsJet)
mapPlot(coastlineWorld, projection="+proj=moll", grid=FALSE)
mapImage(lon, lat, SST, col=oce.colorsJet, zlim=Tlim)
mapPolygon(coastlineWorld, col='gray')