Nothing
# Created by Mao-Gui Hu (humg@lreis.ac.cn), 2013.10.
# 2d EMD for irregular geographical polygons
require(sp)
require(fields)
# Local min/max extremes
LocalExtremes <- function(spXXDf, nbList, boundary, bPlateauAndBeam=TRUE)
{
# spXXDf: SpatialPointsDataFrame or SpatialPolygonsDataFrame
# nbList: neighbours of each objects in spXXDf [list]
# boundary: indicate whether an object in spXXDf is located at boundary [vector]
# bPlateauAndBeam: whether plateau/pan and beam/channel should be treated as extremes
if(!bPlateauAndBeam) {
locMax <- c()
locMin <- c()
for(i in 1:nrow(spXXDf)) {
if(boundary[i]) next
curValue <- as.vector(spXXDf@data$value[i])
nbValue <- as.vector(spXXDf@data$value[nbList[[i]]])
if(all(curValue > nbValue)) locMax <- c(locMax, i)
if(all(curValue < nbValue)) locMin <- c(locMin, i)
}
}
else {
nDiffWithMin <- rep(NA, nrow(spXXDf))
nDiffWithMax <- rep(NA, nrow(spXXDf))
for(i in 1:nrow(spXXDf)) {
if(boundary[i]) next
curV <- as.vector(spXXDf@data$value[i])
minV <- min(c(curV,as.vector(spXXDf@data$value[nbList[[i]]])))
maxV <- max(c(curV,as.vector(spXXDf@data$value[nbList[[i]]])))
nDiffWithMin[i] <- curV - minV
nDiffWithMax[i] <- curV - maxV
}
locMin <- which(nDiffWithMin == 0)
locMax <- which(nDiffWithMax == 0)
}
return(list(locMax=locMax, locMin=locMin))
}
PlotFit <- function(xy, fitrslt, nx=64, ny=64) {
par(mfrow=c(2,2))
quilt.plot(xy, fitrslt$upperEnv, main="Upper envelope", nx=nx, ny=ny)
quilt.plot(xy, fitrslt$lowerEnv, main="Lower envelope", nx=nx, ny=ny)
quilt.plot(xy, fitrslt$meanSurf, main="Mean", nx=nx, ny=ny)
quilt.plot(xy, fitrslt$residual, main="Residual", nx=nx, ny=ny)
}
# Spline interpolation of upper & lower envelopes
SurfaceFit <- function(spPointsDf, locExtr, fmodel="thinplate", fig=FALSE, nx=64, ny=64)
{
# scale coordinates to an appropriate range
ScaleCoord <- function(xy, nmax=100) {
return(xy / max(abs(xy)) * nmax)
}
# cat("Surface fit.\n")
LocalFit <- function()
{
if(fmodel == "cubic" || fmodel == "multiquadric"
|| fmodel == "gaussian" || fmodel == "thinplate") {
xy <- ScaleCoord(spPointsDf@coords)
fitMax <- rbfcreate(xy[locExtr$locMax,], spPointsDf@data$value[locExtr$locMax], rbffun=fmodel)
fitMin <- rbfcreate(xy[locExtr$locMin,], spPointsDf@data$value[locExtr$locMin], rbffun=fmodel)
upperEnv <- as.vector(rbfinterp(fitMax, xy))
lowerEnv <- as.vector(rbfinterp(fitMin, xy))
}
else {
stop("invalid surface fitting method.")
}
meanSurf <- (upperEnv + lowerEnv)/2
residual <- spPointsDf@data$value - meanSurf
fitrslt <- list(#fitMax=fitMax, fitMin=fitMin,
upperEnv=upperEnv, lowerEnv=lowerEnv,
meanSurf=meanSurf, residual=residual)
if(fig) {
PlotFit(spPointsDf@coords, fitrslt, nx=nx, ny=ny)
}
return(fitrslt)
}
tryCatch({
fitrslt <- LocalFit()
return(fitrslt)
},
error=function(x) {
print(x)
return(NULL)
}
)
}
ExtractIMFi2D <- function(spPointsDf, nbList, boundary, tolSift=0.05, minExtrema=5,
fmodel="thinplate", bPlateauAndBeam=TRUE, fig=FALSE, nx=64, ny=64)
{
# spPointsDf: SpatialPointsDataFrame
# nbList: neighbours of each objects in spXXDf [list]
# boundary: indicate whether an object in spXXDf is located at boundary [vector]
# bPlateauAndBeam: whether plateau/pan and beam/channel should be treated as extremes
# cat("Extract IMF.\n")
SD <- Inf
nSift <- 0 # max sift iteration
orgValue <- spPointsDf@data$value
residual <- orgValue
imf <- c()
while(1) {
spPointsDf@data$value <- residual
locExtr <- LocalExtremes(spPointsDf, nbList, boundary, bPlateauAndBeam)
if(length(locExtr$locMax) < minExtrema || length(locExtr$locMin) < minExtrema) break
fit <- SurfaceFit(spPointsDf, locExtr, fmodel=fmodel, fig=fig, nx=nx, ny=ny)
if(is.null(fit) || any(is.nan(fit$residual))) break
SD <- sum((residual-fit$residual)^2)/sum(residual^2)
nSift <- nSift+1
if((SD < tolSift && nSift > 1) || nSift > 5) break
imf <- fit$residual
residual <- imf
}
if(is.null(imf))
return(NULL)
else
residual <- orgValue-imf
if(fig) {
xy <- spPointsDf@coords
par(mfrow=c(1,2))
quilt.plot(xy, imf, main="IMF")
quilt.plot(xy, residual, main="Residual")
}
return(list(imf=imf, residual=residual))
}
MeanIMF <- function(imfs)
{
if(length(imfs) == 1) return(imfs[[1]])
nimf <- rep(0, length(imfs))
for(i in 1:length(imfs)) {
nimf[i] <- NCOL(imfs[[i]]$imf)
}
uniqLen <- sort(unique(nimf))
uniqLenNum <- c()
for(i in uniqLen) {
uniqLenNum <- c(uniqLenNum, sum(nimf==i))
}
selImfs <- imfs[nimf == uniqLen[which.max(uniqLenNum)]]
mnImf <- 0
mnTrend <- 0
for(i in 1:length(selImfs)) {
mnImf <- mnImf+selImfs[[i]]$imf
mnTrend <- mnTrend+selImfs[[i]]$trend
}
mnImf <- mnImf/length(selImfs)
mnTrend <- mnTrend/length(selImfs)
return(list(imf=mnImf, trend=mnTrend))
}
# decompostion of a SpatialPolygonsDataFrame by iBEEMD
# a matrix will be returned:
# the first column is the original value;
# the last column is the decomposed residual (global tredn);
# columns from the second to the last but one is decmompsed IMFs.
iBEEMD <- function(
spPolysDf, # interested spatial data; it is a SpatialPolygonsDataFrame object.
valueField = names(spPolysDf)[1], # interested field value in above spatial data.
nMaxIMF = 10, # maximum number of IMFs (default is 10)
tolSift = 0.05, # sift tolerence (default is 0.05)
neemd = 1000, # EEMD iteration number (default is 1000)
wnsd = 0.05, # standard deviation of added noise; it is a ratio to the standard deviation of above data.
fmodel = "thinplate", # surface fitting function ("thinplate", "gaussian", "cubic", "multiquadric")
fig = TRUE # show decomposed result as an image
)
{
stopifnot(is(spPolysDf, "SpatialPolygonsDataFrame"))
spPolysDfBak <- spPolysDf
minExtrema <- 5 # minimum number of nodes to fit a surface
bPlateauAndBeam <- TRUE # whether plateau/pan and beam/channel should be treated as extremes
spClipBnd <- NULL # if the spPolysDf is a SpatialPointsDataFrame instead of SpatialPolygonsDataFrame,
# a SpatialPolygonsDataFrame dataset will be created by voronoi partition. Then,
# spClipBnd is a SpatialPolygonsDataFrame containing only one polygon to define the boundary.
# reserved for points
bInputPt <- is(spPolysDf, "SpatialPointsDataFrame")
if(bInputPt) {
cat("Create polygons from points ...\n")
spInputPt <- spPolysDf
spPolysDf <- CreateVoronoi(spPolysDf, bDelaunay=FALSE, bClipVor=TRUE, spClipBnd=spClipBnd)
}
spPolysDf@data <- data.frame(id=1:length(spPolysDf), value=spPolysDf@data[,valueField])
noiseSD <- sd(spPolysDf@data$value) * wnsd
if(!bInputPt) {
spPointsDf <- SPolygonsDfToSPointsDf(spPolysDf)
} else {
spPointsDf <- spInputPt
spPointsDf@data <- spPolysDf@data
}
orgLen <- length(spPolysDf)
orgValue <- spPolysDf@data$value
spPolysDf <- createSPComment(spPolysDf)
nbList <- NeighborPolys(spPolysDf)
boundary <- LocatedAtBoundary(spPolysDf, bRemoveHole=TRUE)
cat("EEMD iteration ...\n")
if(noiseSD == 0 || neemd <= 0) neemd <- 1
imfs <- list()
bakValue <- spPointsDf@data$value
pb <- txtProgressBar(min=0, max=neemd, width=50, style=3)
for(k in 1:neemd)
{
# extract imfs
imf <- c()
trend <- c()
spPointsDf@data$value <- bakValue+rnorm(length(spPointsDf),sd=noiseSD)
for(i in 1:nMaxIMF)
{
extimf <- ExtractIMFi2D(spPointsDf, nbList, boundary, tolSift=tolSift, minExtrema=minExtrema,
bPlateauAndBeam=bPlateauAndBeam, fig=FALSE, fmodel=fmodel)
if(is.null(extimf)) break
imf <- cbind(imf, as.vector(extimf$imf))
trend <- extimf$residual
spPointsDf@data$value <- trend
if(sum(trend) == 0) break
}
if(is.null(imf)) next
imfs[[length(imfs)+1]] <- list(imf=imf, trend=as.vector(trend))
setTxtProgressBar(pb, k)
}
close(pb)
if(length(imfs) == 0) {
warning("no imf decomposed.")
return(NULL)
}
if(!is.null(imfs))
mnImfTrend <- MeanIMF(imfs)
imfDf <- as.data.frame(mnImfTrend$imf)
colnames(imfDf) <- paste("imf", 1:ncol(imfDf), sep="")
spPolysDfBak@data <- data.frame(original=orgValue, imfDf, trend=mnImfTrend$trend)
names(spPolysDfBak@data)[1] <- valueField
if(fig) {
cat("Plotting ...\n")
print(spplot(spPolysDfBak, col.regions = terrain.colors(256)))
}
return(spPolysDfBak)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.