[R-sig-Geo] Madogram, rodogram, semi-variogram of satellite imagery

[Ashton Shortridge]

Hi Chethan,

The fields library does variograms efficiently on images. Check out the 
vgram.matrix function in there.

Cressie's robust variogram is already implemented, while madograms (or 
rodograms) are not implemented. However, the vgram.matrix function is basic R. 
I modified it readily to calculate madograms (R file with implemented function 
is attached); you could implement rodograms as easily. As far as I know cross-
variograms are not implemented in fields.

Here is an example using the modified function (attached):
library(fields)
data(lennon)
image(lennon)

source("madogram_fields.R")   # loads a modified function called vgram2.matrix()
out2 <- vgram2.matrix(lennon)
plot(out2$d.full, out2$madgram, xlab="separation distance", ylab="madogram")
title(main="Madogram for Lennon Image")

Hope this helps,

Ashton


On 2011-01-22, Chethan S, wrote:
> Greetings!
> 
> I am cross-posting this query here as suggested by R-Help list.
> 
> Is there in any package for R which can help me generate madogram,
> rodogram, semi-variogram, cross variogram from landsat imagery. I intend
> to select portions of large imagery (i.e., obtain subsets) and generate
> the above said texture layers. A simple Google search led to me a result
> suggesting SpatialExtremes for generating madograms. However I am not sure
> about its capability to generate one for image layers. I am a beginner
> level user in R and that's why I am posting this question to give me a
> quick start.
> 
> Thanks and regards,
> 
> Chethan S.


-----
Ashton Shortridge
Associate Professor			ashton at msu.edu
Dept of Geography			http://www.msu.edu/~ashton
235 Geography Building		ph (517) 432-3561
Michigan State University		fx (517) 432-1671
-------------- next part --------------

vgram2.matrix <- function (dat, R = 5, dx = 1, dy = 1) 
{
    # This function is a slight modification of vgram.matrix() in the 
    # fields library to calculate a madogram on matrix (image) data.
    # In fact, only two lines were added! Check help on vgram.matrix 
    # for more information. Requires the fields library. 
    # Written by Ashton Shortridge, 1/2011
    SI <- function(ntemp, delta) {
        n1 <- 1:ntemp
        n2 <- n1 + delta
        good <- (n2 >= 1) & (n2 <= ntemp)
        cbind(n1[good], n2[good])
    }
    N <- ncol(dat)
    M <- nrow(dat)
    m <- min(c(round(R/dx), M))
    n <- min(c(round(R/dy), N))
    ind <- rbind(as.matrix(expand.grid(0, 1:n)), as.matrix(expand.grid(1:m, 
        0)), as.matrix(expand.grid(c(-(m:1), 1:m), 1:n)))
    d <- sqrt((dx * ind[, 1])^2 + (dy * ind[, 2])^2)
    good <- (d > 0) & (d <= R)
    ind <- ind[good, ]
    d <- d[good]
    ind <- ind[order(d), ]
    d <- sort(d)
    nbin <- nrow(ind)
    holdVG <- rep(NA, nbin)
    holdMG <- rep(NA, nbin)
    holdRVG <- rep(NA, nbin)
    holdN <- rep(NA, nbin)
    for (k in 1:nbin) {
        MM <- SI(M, ind[k, 1])
        NN <- SI(N, ind[k, 2])
        holdN[k] <- length(MM) * length(NN)
        BigDiff <- (dat[MM[, 1], NN[, 1]] - dat[MM[, 2], NN[, 
            2]])
        holdVG[k] <- mean(0.5 * (BigDiff)^2)   # variogram
        holdMG[k] <- mean(0.5 * abs(BigDiff))  # madogram
        holdRVG[k] <- mean(abs(BigDiff)^0.5)   # Cressie & Hawkins
    }
    holdRVG <- 0.5 * (holdRVG^4)/(0.457 + 0.494 * holdN)
    top <- tapply(holdVG * holdN, d, FUN = "sum")
    bottom <- tapply(holdN, d, FUN = "sum")
    dcollapsed <- as.numeric(names(bottom))
    vgram <- top/bottom
    dimnames(vgram) <- NULL
    list(vgram = vgram, d = dcollapsed, ind = ind, d.full = d, 
        vgram.full = holdVG, robust.vgram = holdRVG, madgram = holdMG, 
        N = holdN, dx = dx, dy = dy)
}