Re: [R-sig-Geo] raster: stackApply problems..
Thank you Jon!
In fact, that's how I thought it worked.
And that's how it worked for me all the time!
But recently, doing some manual checks on some indices I couldn't
confirm it ...
I tried to replicate the problem and my workflow with test data
(https://gist.github.com/kokkytos/5d554b5a725bb48d2189e2d1fa0e2206) but
stackApply works fine!
However, in my real data/workflow when I try the same procedure, I have
issues with the returned names of stackApply (indexes of names do not
match).
That's the result of the real data to see what I mean, name indices
doesn't match:
> ver_median (my verification with alternative way)
class : RasterStack
dimensions : 300, 300, 9, 7 (nrow, ncol, ncell, nlayers)
resolution : 500, 500 (x, y)
extent : 379354, 529354, 4132181, 4282181 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=50 +y_0=0
+ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs
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 : 297.5, 311.0, 313.0, 468.0, 290.0, 302.0, 307.0
> stackapply_median (stackapply calculations)
class : RasterBrick
dimensions : 300, 300, 9, 7 (nrow, ncol, ncell, nlayers)
resolution : 500, 500 (x, y)
extent : 379354, 529354, 4132181, 4282181 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=50 +y_0=0
+ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs
source : /home/leonidas/tmp/r_tmp_2019-11-22_185242_4475_86751.grd
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
I'll look again at my code...
I apologize for the inconvenience.
On 11/22/19 6:31 PM, jon.sko...@ec.europa.eu wrote:
>
> 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 t
Re: [R-sig-Geo] raster: stackApply problems..
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] Comparing distance among point pattern events
Hi,
Great question, and clear example.
The first problem:
ACd<-pairdist(A) instead of ACd <- pairdist(AC)
BUT
pairdist() is the wrong function: that calculates the mean distance
between ALL points, A to A and C to C as well as A to C.
You need crossdist() instead.
The most flexible approach is to roll your own permutation test. That
will work even if B and C are different sizes, etc. If you specify the
problem more exactly, there are probably parametric tests, but I like
permutation tests.
library(spatstat)
set.seed(2019)
A <- rpoispp(100) ## First event
B <- rpoispp(50) ## Second event
C <- rpoispp(50) ## Third event
plot(A, pch=16)
plot(B, col="red", add=T)
plot(C, col="blue", add=T)
ABd<-crossdist(A, B)
ACd<-crossdist(A, C)
mean(ABd)
# 0.5168865
mean(ACd)
# 0.5070118
# test the hypothesis that ABd is equal to ACd
nperm <- 999
permout <- data.frame(ABd = rep(NA, nperm), ACd = rep(NA, nperm))
# create framework for a random assignment of B and C to the existing points
BC <- superimpose(B, C)
B.len <- npoints(B)
C.len <- npoints(C)
B.sampvect <- c(rep(TRUE, B.len), rep(FALSE, C.len))
set.seed(2019)
for(i in seq_len(nperm)) {
B.sampvect <- sample(B.sampvect)
B.perm <- BC[B.sampvect]
C.perm <- BC[!B.sampvect]
permout[i, ] <- c(mean(crossdist(A, B.perm)), mean(crossdist(A, C.perm)))
}
boxplot(permout$ABd - permout$ACd)
points(1, mean(ABd) - mean(ACd), col="red")
table(abs(mean(ABd) - mean(ACd)) >= abs(permout$ABd - permout$ACd))
# FALSE TRUE
# 573 426
sum(abs(mean(ABd) - mean(ACd)) >= abs(permout$ABd - permout$ACd)) / nperm
# 0.4264264
The difference between ACd and ABd is indistinguishable from that
obtained by a random resampling of B and C.
Sarah
On Fri, Nov 22, 2019 at 8:26 AM ASANTOS via R-sig-Geo
wrote:
>
> Dear R-Sig-Geo Members,
>
> I have the hypothetical point process situation:
>
> library(spatstat)
> set.seed(2019)
> A <- rpoispp(100) ## First event
> B <- rpoispp(50) ## Second event
> C <- rpoispp(50) ## Third event
> plot(A, pch=16)
> plot(B, col="red", add=T)
> plot(C, col="blue", add=T)
>
> I've like to know an adequate spatial approach for comparing if on
> average the event B or C is more close to A. For this, I try to make:
>
> AB<-superimpose(A,B)
> ABd<-pairdist(AB)
> AC<-superimpose(A,C)
> ACd<-pairdist(A)
> mean(ABd)
> #[1] 0.5112954
> mean(ACd)
> #[1] 0.5035042
>
> With this naive approach, I concluded that event C is more close of A
> that B. This sounds enough for a final conclusion or more robust
> analysis is possible?
>
> Thanks in advance,
>
> Alexandre
>
--
Sarah Goslee (she/her)
http://www.numberwright.com
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Comparing distance among point pattern events
Dear R-Sig-Geo Members, I have the hypothetical point process situation: library(spatstat) set.seed(2019) A <- rpoispp(100) ## First event B <- rpoispp(50) ## Second event C <- rpoispp(50) ## Third event plot(A, pch=16) plot(B, col="red", add=T) plot(C, col="blue", add=T) I've like to know an adequate spatial approach for comparing if on average the event B or C is more close to A. For this, I try to make: AB<-superimpose(A,B) ABd<-pairdist(AB) AC<-superimpose(A,C) ACd<-pairdist(A) mean(ABd) #[1] 0.5112954 mean(ACd) #[1] 0.5035042 With this naive approach, I concluded that event C is more close of A that B. This sounds enough for a final conclusion or more robust analysis is possible? Thanks in advance, Alexandre -- Alexandre dos Santos Geotechnologies and Spatial Statistics applied to Forest Entomology Instituto Federal de Mato Grosso (IFMT) - Campus Caceres Caixa Postal 244 (PO Box) Avenida dos Ramires, s/n - Distrito Industrial Caceres - MT - CEP 78.200-000 (ZIP code) Phone: (+55) 65 99686-6970 / (+55) 65 3221-2674 Lattes CV: http://lattes.cnpq.br/1360403201088680 OrcID: orcid.org/-0001-8232-6722 ResearchGate: www.researchgate.net/profile/Alexandre_Santos10 Publons: https://publons.com/researcher/3085587/alexandre-dos-santos/ -- ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
