Re: [R-sig-Geo] raster: stackApply problems..

[Jon.SKOIEN]

Leonidas,

I see that you are not happy with the output, but it is not so clear what you 
actually expect to see.


If you use stackApply directly, the indices are used in the names. Layer 1 and 
8 belong to the group with index 4. It is the first group in the list of 
indexes, so the first layer of the output is then referred to as index_4. Then 
comes index_5 with layers 2, 10 and 15 of your input. The order of these names 
will follow the order of the first appearance of your indices. The indices gets 
lost with the use of clusterR, so it gives you the same output, but with names 
layer.1 - layer.7.


You could change the names of the result from clusterR with:

names(ResClusterR) = paste0("index_", unique(indices))


If you want your result (from stackApply or clusterR) to follow the order of 
your indices, you should be able to get this with:


sResClusterR = ResClusterR[[order(names(ResClusterR))]]


Does this help you further?

Jon




--
Jon Olav Skøien
European Commission
Joint Research Centre – JRC.E.1
Disaster Risk Management Unit
Building 26b 1/144 | Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 267
jon.sko...@ec.europa.eu
 Tel:  +39 0332 789205 Disclaimer: Views expressed in this email are those of 
the individual and do not necessarily represent official views of the European 
Commission.




From: R-sig-Geo  on behalf of Leonidas Liakos 
via R-sig-Geo 
Sent: 21 November 2019 08:52
To: Ben Tupper; r-sig-geo@r-project.org
Subject: Re: [R-sig-Geo] raster: stackApply problems..

Unfortunately the names are not always in ascending order. This is the
result of my data.

names  : index_4, index_5, index_6, index_7, index_1, index_2, index_3
min values :   3,   3,   3,   3,   3,   3,   3
max values :   307.0,   297.5,   311.0,   313.0,   468.0,   290.0,   302.0

And worst of all, it is not a proper match with indices.

If I run it with clusterR then the result is different:

names  : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7
min values :   3,   3,   3,   3,   3,   3,   3
max values :   307.0,   297.5,   311.0,   313.0,   468.0,   290.0,   302.0


The solution is to reorder the layers of the stack so that the
stackApply indices are in ascending order e.g. 1,1,1,2,2,2,3,3,3 ...

My indices of my data was like that:

4 5 6 7 1 2 3 4 5 6 7 1 2 3 5 6 7

I've reported this behavior here
https://urldefense.com/v3/__https://github.com/rspatial/raster/issues/82__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G-KjNZJC$


On 11/20/19 3:05 PM, Ben Tupper wrote:
> Hi,
>
> That is certainly is unexpected to have two different naming styles.
> It's not really solution to take to the bank, but you could simply
> compose your own names assuming that the layer orders are always
> returned in ascending index order.
> Would that work for you
>
> ### start
> library(raster)
>
> # Compute layer names for stackApply output
> #
> # @param index numeric, 1-based layer indices used for stackApply function
> # @param prefix character, prefix for names
> # @return character layers names in index order
> layer_names <- function(index = c(2,2,3,3,1,1), prefix = c("layer.",
> "index_")[1]){
>   paste0(prefix, sort(unique(index)))
> }
>
> indices <- c(2,2,3,3,1,1)
>
> r <- raster()
> values(r) <- 1
> # simple sequential stack from 1 to 6 in all cells
> s <- stack(r, r*2, r*3, r*4, r*5, r*6)
> s
>
> beginCluster(2)
> res <- clusterR(s, stackApply, args = list(indices=indices, fun = mean))
> raster::endCluster()
> names(res) <- layer_names(indices, prefix = "foobar.")
> res
>
> res2 <- stackApply(s, indices, mean)
> names(res2) <- layer_names(indices, prefix = "foobar.")
> res2
> ### end
>
>
> On Wed, Nov 20, 2019 at 1:36 AM Leonidas Liakos via R-sig-Geo
>  wrote:
>> This is not a reasonable solution. It is not efficient to run stackapply
>> twice to get the right names. Each execution can take hours.
>>
>>
>> Στις 20/11/2019 3:30 π.μ., ο Frederico Faleiro έγραψε:
>>> Hi Leonidas,
>>>
>>> both results are in the same order, but the name is different.
>>> You can rename the first as in the second:
>>> names(res) <- names(res2)
>>>
>>> I provided an example to help you understand the logic.
>>>
>>> library(raster)
>>> beginCluster(2)
>>> r <- raster()
>>> values(r) <- 1
>>> # simple sequential stack from 1 to 6 in all cells
>>> s <- stack(r, r*2, r*3, r*4, r*5, r*6)
>>> s
>>> res <- clusterR(s, stackApply, args = list(indices=c(2,2,3,3,1,1), fun =
>>> mean))
>>> res
>>> res2 <- stackApply(s, c(2,2,3,3,1,1), mean)
>>> res2
>>> dif <- res - res2
>>> # exatly the same order because the difference is zero for all layers
>>> dif
>>> # rename
>>> names(res) <- names(res2)
>>>
>>> Best regards,
>>>
>>> Frederico Faleiro
>>>
>>> On Tue, Nov 19, 2019 at 

Re: [R-sig-Geo] Warning in fit.variogram: value out of range in 'bessel k' | automap and gstat packages

[Jon.SKOIEN]

Frederico,


I don't think you need to worry about this warning in this case. 
autofitVariogram tests a lot of variogram models (and different kappa-values) 
in the search for the best one, and then fit.variogram tries to optimize the 
parameters based on these. It seems the warnings were generated in the C-code 
of gstat, based on some of the tested kappa-values for the Matern variogram 
(Stein implementation). It also seems that these kappa values did not give good 
fits to the sample variogram, so you did not get the warning from the best fit 
variogram in this case.


Hope this helps,

Jon



--
Jon Olav Sk�ien
European Commission
Joint Research Centre � JRC.E.1
Disaster Risk Management Unit
Building 26b 1/144 | Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 267
jon.sko...@ec.europa.eu
 Tel:  +39 0332 789205 Disclaimer: Views expressed in this email are those of 
the individual and do not necessarily represent official views of the European 
Commission.


From: R-sig-Geo  on behalf of Frederico 
Faleiro 
Sent: 09 September 2019 21:28:32
To: RsigGeo
Subject: [R-sig-Geo] Warning in fit.variogram: value out of range in 'bessel k' 
| automap and gstat packages

Dear list,

I am trying make an interpolation with the function autoKrige from the
package automap, but I received more than 50 times the following warning
message:

In fit.variogram(experimental variogram, model = vgm(psill = psill, ... :
value out of range in 'bessel k'

I have tryed identify it in the code of fit.variogram and vgm from the
gstat package, but they do not have any warning about it in the code. Then
I imagine the warning is from another function.
I fitted the variogram automatically in the automap and after I fitted in
the gstat to check if there is the same warning for that parameters, but in
the last I do not have any warning. When I check the variogram model, it
has a good fit apparently. I provided the reproductible example below.

Do you know the consequences of this warning and how can I avoid it? I
would like to use the automap package because I need interpolate many other
files like this.

# example
library(automap)
library(sp)

# download the data in:
https://urldefense.proofpoint.com/v2/url?u=https-3A__drive.google.com_file_d_1L33-5FQgpNMdwvWff-5FuMrwhl-5FKCajrtlVe_view-3Fusp-3Dsharing=DwICAg=8NwulVB6ucrjuSGiwL_ckQ=DA6jV5X00Z1YpyMh4OS79xeSQGErBqY83CBz841DNnU=YlpT3yKlbNVLvXMp08O8X1LkniHGXTzJeLDHLTotwOI=U2z98vHb1mpKhKJ9UXA0RUrPi6r6IXibq1cuuyB0jWY=
# read the data
pr <- read.table('pr_monC.txt', header = TRUE, sep = "\t")
coordinates(pr) <- ~long+lat
# read the grid
gr <- read.table('grid.txt' , header = TRUE, sep = '\t')
gridded(gr) <-  ~long+lat

# automatic fit variogram with automap
av <- autofitVariogram(jan~1, pr)
# warnings message:   In fit.variogram(experimental variogram, model =
vgm(psill = psill, ... : value out of range in 'bessel k'
plot(av) # get in the graph the parameters to fit the variogram

# fit variogram in gstat with the model parameters from autofitVariogram to
check if there is any warning message too
ve <- variogram(jan~1, pr)
v <- fit.variogram(ve, vgm(psill = 9723, model = "Ste", range = 20, nugget
= 194, kappa = 0.8))
# it does not produce any warning
plot(ve, v)

Best regards,

Frederico

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dsig-2Dgeo=DwICAg=8NwulVB6ucrjuSGiwL_ckQ=DA6jV5X00Z1YpyMh4OS79xeSQGErBqY83CBz841DNnU=YlpT3yKlbNVLvXMp08O8X1LkniHGXTzJeLDHLTotwOI=wzYnAjR9sQtl23JiTMc5f1m79m8uQgUt9DSCj3wuTeA=

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Problems with gstat gaussian variogram and cross validation when fix nugget at 0

[Jon.SKOIEN]

Stefano,

Fixing the nugget at zero for the Gaussian varigoram can often cause numerical 
problems in kriging. If you have two observations that are close in space, 
their modelled semivariance will be so small that you get numerical issues with 
the covariance matrix. In some cases it will be singular, in your case it could 
be inverted, but the absolute values of the weights will be very high for some 
of the cross-validation locations. It seems you have two stations that are very 
close to each other, if you remove one of these, you can see how the results 
will be more similar (still not the same) to the model with a small nugget 
effect.

Best,
Jon

--
Jon Olav Skøien

European Commission
Joint Research Centre – JRC.E.1
Disaster Risk Management Unit
Building 26b 1/144 | Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 267
jon.sko...@ec.europa.eu
Tel:  +39 0332 789205
Disclaimer: Views expressed in this email are those of the individual and do 
not necessarily represent official views of the European Commission.



From: R-sig-Geo [r-sig-geo-boun...@r-project.org] on behalf of Stefano 
Menichetti [s.meniche...@arpat.toscana.it]
Sent: 12 October 2018 16:36
To: r-sig-geo@r-project.org
Subject: [R-sig-Geo] Problems with gstat gaussian variogram and cross 
validation when fix nugget at 0

Dear all, I have a question about fixing nugget at 0 in a Gaussian
variogram using *gstat *package

I have encountered a serious problem on a variogram of piezometric data
that give "strange" results if I want to fix nugget to 0.

This is the file and how to import:

/library(gstat)/
//file =>
https://drive.google.com/open?id=1AOilEhR8G8YANVuFIBUcmEg-6UG1Ajrr//

/lay1 =� read.delim(paste(path,*file*,sep=""),sep="\t",dec=",")//
//SpData = na.exclude(subset(lay1, select = c("X","Y","Observed"))) //
//colnames(SpData)=c("X","Y","Z")/

The initial variogram:

/coordinates(SpData) = ~X+Y//
//var = variogram(Z~1,SpData)/

show a slow initial rise and so suggest a *Gaussian *type.

After a first variogram with:

/m = vgm(200,"Gau",3000,0)//
//
//plot(var, model=m, //
//main = paste(m$model[2],"//
//Nug=",m$psill[1],", Sill=",round(m$psill[2]),",
Range=",round(m$range[2]),sep ="")//
//)/

I have obtained a fitted variogram like this,

/m = fit.variogram(var, m)/

with a reasonable cross-validation results:

/krg_cv = krige.cv(Z~1, SpData, model = m)//
//plot(krg_cv$observed,krg_cv$var1.pred, xlab="Z observed", ylab="Z
predicted")//
//abline(a=0,b=1)//
//paste("media residui =",mean(krg_cv$residual))//
//paste("RMS residui =",mean(krg_cv$residual^2)^0.5)/

Resulting map it is very smooth, despite the small, it seems, nugget and
so *I have tried to fix nugget at 0* (exact interpolation).

But cross validation results with

/m = vgm(216,"Gau",3640,*0*)/

are dramatically wrong and also very different from cv results with

/m = vgm(216,"Gau",3640,*1e-3*)/

What happens ? I wrong in something ? Why nugget paraneter at 0 or at
small value influnce so hard the results ?

Thank you very much in advance,


Stefano


Dott. Geol. Stefano Menichetti
=
ARPAT - Agenzia per la Protezione dell'Ambiente della Toscana
Settore Tecnico SIRA - Sistema Informativo Regionale Ambientale via Nicola 
Porpora 22 - 50144 Firenze
tel. 0553206333 cell. 3668217978 (3383550147) fax 0553206410
s.meniche...@arpat.toscana.it  http:\\www.arpat.toscana.it  
http:\\sira.arpat.toscana.it


[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo