Creating simple masks

One of the required components for a tripEstimation model is a "mask". This is implemented as a lookup function which provides a test for coordinates in a matrix - "which of these coordinates are allowed by the mask". The mask may be as simple as a bounding box (i.e. no other data is required apart from four values xmin, xmax, ymin, ymax) or as complex as a spatio-temporal mask comparing temperature ranges from a tag to satellite-measured ocean temperatures.

The following examples show some ways of building the lookup functions that are required by the tripEstimation models.

Bounding box lookup mask

The only inputs required in this simple case are the four values describing the min and max coordinates.

lon.min <- 100
lon.max <- 170
lat.min <- -45
lat.max <- 80

The following function returns a lookup function that does a test on whether coordinates fall within these bounds.

## lookup function that encapsulates the bounding box definition, this function
## takes a matrix of coordinates of [lon,lat] and returns a logical
## vector for each lon,lat pair of whether it is allowed by the mask

makelookup <- function(xlim, ylim) {
    function(x) {
        x[,1] >= xlim[1] & x[,1] <= xlim[2] & x[,2] >= ylim[1] & x[,2] <= ylim[2]
    }
}

We run makelookup with the values defined above, and it returns a function that can be used directly in the model templates (e.g. solar.model()) in tripEstimation.

lookup <- makelookup(xlim = c(lon.min, lon.max), ylim = c(lat.min, lat.max))

Raster mask from world land polygons

Define a land mask where the lookup function tests whether a coordinate falls on the land area versus the sea.

## raster mask from built-in polygons

library(raster)
library(maptools)

## load in the world polygons, a built-in data set in maptools
data(wrld_simpl)

## build a raster with a given dimension
r <- raster(nrows = 250, ncols = 250, xmn = lon.min, xmx = lon.max,
            ymn = lat.min, ymx = lat.max, crs = proj4string(wrld_simpl))

## populate the raster with values from the polygons
r <- rasterize(wrld_simpl, r)

## process the values in order to be a mask
r <- is.na(r)  ## sea is now 1

## degrade this raster to an R list with components x, y, z which is
## what tripEstimation expects (tripEstimation could be trained to do this)
mask <- as.image.SpatialGridDataFrame(as(r, "SpatialGridDataFrame"))

## land and not the ocean :)
mask$z <- !mask$z

## now create the lookup function
library(tripEstimation)
lookup <- mkLookup(mask, by.segment = FALSE)

## convinced?

##image(mask, main = "click on the map")
##while(TRUE) {
##    xy <- locator(1)
##    text(xy$x, xy$y, lab = lookup(cbind(xy$x, xy$y)))
##}

Reading in data from shapefiles, databases, other polygon formats

TODO

Create a mask directly from a topographic raster data source

For this example we use the Etopo1 work data set of topography.

The data is available from here: http://www.ngdc.noaa.gov/mgg/global/

The specific file below used is downloaded from here and unzipped: http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/cell_registered/georeferenced_tiff/ETOPO1_Ice_c_geotiff.zip

## raster mask from ETOPO1
library(raster)
library(rgdal)

r <- raster("ETOPO1_Ice_g_geotiff.tif")


## crop the raster
r <- crop(r, extent(lon.min, lon.max, lat.min, lat.max))

## process the logic of the mask
r <- r < 500 & r > -2000

## degrade this raster to an R list with components x, y, z which is
## what tripEstimation expects (tripEstimation could be trained to do this)
mask <- as.image.SpatialGridDataFrame(as(r, "SpatialGridDataFrame"))

## if you want the land and not the ocean :)
##mask$z <- !mask$z

## now create the lookup function
library(tripEstimation)
lookup <- mkLookup(mask, by.segment = FALSE)