[R-sig-Geo] Grid projections and spatial matrices

[Colin Beale]

Hi,

Two questions really (elaborated below): (1) Has anyone created what is
essentially a spherical (or at least hemispherical) matrix structure
that I can use normal R functions like apply() on? and (2) has anyone
coded the procedure described in the ref below to calculate effective
sample sizes of correlations between two spatial matrices? As the first
author is clearly an R user I had expected to find something, but quite
a lot of searching has so far drawn a blank.

Clifford, P., Richardson, S. & Hemon, D. (1989). Assessing the
significance of the correlation between two spatial processes.
Biometrics, 45, 123-134.

As the first question is the most involved, this is the background: I'm
looking at presence / absence data for a range of species. Up to now, I
have mainly worked with distributions easily mapable on a regular grid.
I've developed a number of functions that use arrays as a way of
representing these data: y & x references in the first two dimensions,
species ID in a third. Several of these functions compute turnover in
species distirbution - they check neighbouring boxes within the grid for
presence / absence of the same species and compute statistics based on
this. This is all very well when the basic projection is a simple 2d
matrix. Recently, however, I have started using larger scale datasets
where appropriate projections need to be used - most of this data is
projected onto a 10x10km UTM based grid reference system (in practice,
the Military Grid Reference System). I need to adapt my functions to
recognise a dataset that has this projection - ie, that cells in
different sectors can be neighbours. I appreciate that this might
require a complete rewriting of my functions and a new way of thinking
about these data, but the most obvious solution to me is simply to
create a matrix or array type object that is sort-of "zipped" together
along the edges of the sector, to retain appropriate neighbour
identification. As an example of the sort of basic calculation involved,
I include below a small function below that calcualtes one parameter of
interest - related to the number of species that a species if present in
the focal quadrat, an in its neighbouring ones that I currently use on
my flat 2d array type data - I've also included some toy distributions
so you can see how it works. Any suggestions/pointers would be very
welcome (and I'm sure the original code could be more elegant too!).

Thanks very much,

Colin

# toy dataset
set.seed (1)
distMap <- array (sample (c (0, 0, 1), 120, rep = T), dim = c (6, 5,
4))

##example function - require package(magic) to run.

localTurnover <- function (spmap = distMap) {  
  mapx    <- dim (spmap)[2]                             # extracts map
dimensions and start 
  mapy    <- dim (spmap)[1]                             # up
parameters
  nsp     <- dim (spmap)[3]
  library (magic)                
  temp <- apad(apad(spmap, c(1,1,0)), c(1,1,0), post = FALSE)  # pads
array
  temp[,c (1, dim(temp)[2]),] <- NA                            # makes
additional rows/cols = NA
  temp[c (1, dim(temp)[1]),,] <- NA                            
  spmap[!is.na(spmap)] <- 1                                    # makes
all values in spmap = 1
  turnoverMap <- spmap

  locala <- rep (1, length = 8)
  localb <- rep (1, length = 8)
  localc <- rep (1, length = 8)
 
  for (NSP in 1:nsp) {
    for (iy in 1:mapy) {
      for (ix in 1:mapx) {
        if (!is.na (temp[iy+1,ix+1,NSP])) {
          if (temp[iy+1,ix+1,NSP] == 0) turnoverMap[iy, ix, NSP] <- 0
else {
            if (is.na (temp[iy    , ix    , NSP])) locala[1] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy    ,ix     , NSP]
== 0) locala[1] <- 0
            }
            if (is.na (temp[iy + 1, ix    , NSP])) locala[2] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy + 1,ix     , NSP]
== 0) locala[2] <- 0
            }
            if (is.na (temp[iy + 2, ix    , NSP])) locala[3] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy + 2,ix     , NSP]
== 0) locala[3] <- 0
            }
            if (is.na (temp[iy    , ix + 1, NSP])) locala[4] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy    , ix + 1, NSP]
== 0) locala[4] <- 0
            }   
            if (is.na (temp[iy    , ix + 2, NSP])) locala[5] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy    , ix + 2, NSP]
== 0) locala[5] <- 0
            }
            if (is.na (temp[iy + 1, ix + 2, NSP])) locala[6] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy + 1, ix + 2, NSP]
== 0) locala[6] <- 0
            }
            if (is.na (temp[iy + 2, ix + 1, NSP])) locala[7] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy + 2, ix + 1, NSP]
== 0) locala[7] <- 0
            }
            if (is.na (temp[iy + 2, ix + 2, NSP])) locala[8] <- NA else
{
              if (temp[iy+1,ix+1,NSP] == 1 & temp[iy + 2, ix + 2, NSP]
== 0) locala[8] <- 0
            }
            turnoverMap[iy,ix,NSP] <- mean (locala, na.rm = TRUE)
            locala <- rep (1, length = 8)
          }
        }  
      }
    }
  }
  return (apply(turnoverMap, c(1,2), sum))
}

localTurnover()

Dr. Colin Beale
Spatial Ecologist
The Macaulay Institute
Craigiebuckler
Aberdeen
AB15 8QH
UK

Tel: 01224 498200 ext. 2363
Fax: 01224 311556
Email: c.beale at macaulay.ac.uk