knitr::opts_chunk$set(fig.path = "figure/newmnfTestima_")
require(spacetime)
require(raster)
require(rgdal)
dirdata <- "/media/alobo/LACIE500/MNFtests/data/midata"
ori <- brick(file.path(dirdata,"mini.img"))
noise <- brick(file.path(dirdata,"noise.envi"))
noised <- brick(file.path(dirdata,"mininoise.envi"))
noisecov3 <- layerStats(noise,'cov')$covariance
#it would be good having the maf method also to estimate noise
#over an specific part of the input data. i.e.
#noisecov3 <- maf(layerStats(crop(noised,polygon)),'cov')$covariance)
#were polygon would define a suitable region
testmnf <- mnf(as(noised, 'SpatialGridDataFrame'),Sigma.Noise=noisecov3)
str(testmnf)
## List of 6
## $ values : num [1:6] 0.9128 0.7892 0.572 0.1556 0.0199 ...
## $ sdev : num [1:6] 0.955 0.888 0.756 0.394 0.141 ...
## $ rotation: num [1:6, 1:6] -0.4695 0.81 -0.3429 -0.0474 0.0523 ...
## $ x :Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
## .. ..@ data :'data.frame': 383130 obs. of 6 variables:
## .. .. ..$ MNF1: num [1:383130] 0.00827 0.00196 -0.00983 -0.00275 -0.00632 ...
## .. .. ..$ MNF2: num [1:383130] 0.0295 0.0283 0.0294 0.0129 0.0143 ...
## .. .. ..$ MNF3: num [1:383130] 0.0525 0.0475 0.0614 0.0379 0.0298 ...
## .. .. ..$ MNF4: num [1:383130] 0.0425 0.0593 0.0387 0.05 0.0553 ...
## .. .. ..$ MNF5: num [1:383130] -0.0409 -0.0583 -0.0375 -0.0379 -0.04 ...
## .. .. ..$ MNF6: num [1:383130] 0.435 0.441 0.466 0.432 0.422 ...
## .. ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
## .. .. .. ..@ cellcentre.offset: Named num [1:2] 299452 4622122
## .. .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
## .. .. .. ..@ cellsize : num [1:2] 15 15
## .. .. .. ..@ cells.dim : int [1:2] 774 495
## .. ..@ bbox : num [1:2, 1:2] 299445 4622115 311055 4629540
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "s1" "s2"
## .. .. .. ..$ : chr [1:2] "min" "max"
## .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
## .. .. .. ..@ projargs: chr "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
## $ scale : logi FALSE
## $ center : logi FALSE
## - attr(*, "class")= chr [1:2] "mnf" "prcomp"
plot(testmnf$values-1)
testmnfx <- testmnf$x
testmnfx <- brick(testmnfx)
#writeRaster(testmnfx,file="testmnfx",format="GTiff")
Display
op <-par()
par(mar=c(.1,.1,1,.1),cex.lab=.95,cex.axis=.9,mgp=c(2,.7,0),tcl=-.3)
plot(stretch(testmnfx,minq=0.1, maxq=0.9),col=grey.colors(64),
legend=FALSE,axes=FALSE)
Note they are not ordered as the eigenvalues Comparing to http://diapiro.ictja.csic.es/alobo/blog/CompareOrig2invMNF_log.html) I would say they are 5,6,4,3,2,1 Also, I think that svd() is better than eigen (see help page of princomp) in MNF.
I reorder to apply the inverse
newx <- testmnfx[[c(5,6,4,3,2,1)]]
newx[[6]] <- cellStats(newx[[6]],'mean')
newx[[5]] <- cellStats(newx[[5]],'mean')
newx <- calc(newx,function(x)x%*%solve(testmnf$rotation[,c(5,6,4,3,2,1)]))
#writeRaster(newx,file="newx",format="GTiff")
plot(stretch(newx,minq=0.1, maxq=0.9),col=grey.colors(64),
legend=FALSE,axes=FALSE)
TBC… I will also improve the image display