Re: raster: stackApply problems..

classic Classic list List threaded Threaded
17 messages Options
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
I run the example with clusterR:

no_cores <- parallel::detectCores() -1
raster::beginCluster(no_cores)
?????? res <- raster::clusterR(inp, raster::stackApply, args =
list(indices=c(2,2,3,3,1,1),fun = mean))
raster::endCluster()

And the result is:

> res
class?????????? : RasterBrick
dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
resolution : 1, 1?? (x, y)
extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source???????? : memory
names?????????? : layer.1, layer.2, layer.3
min values :???????? 1.5,???????? 3.5,???????? 5.5
max values :???????? 1.5,???????? 3.5,???????? 5.5??


layer.1, layer.2, layer.3 (?)

So what corrensponds to what?


If I run:

res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)

The result is:

> res2
class      : RasterBrick
dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : index_2, index_3, index_1
min values :     1.5,     3.5,     5.5
max values :     1.5,     3.5,     5.5

There is no consistency with the names of the output and obscure correspondence with the indices in the case of clusterR


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

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 4:15 PM Leonidas Liakos via R-sig-Geo <
[hidden email]> wrote:

> I run the example with clusterR:
>
> no_cores <- parallel::detectCores() -1
> raster::beginCluster(no_cores)
> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
> list(indices=c(2,2,3,3,1,1),fun = mean))
> raster::endCluster()
>
> And the result is:
>
> > res
> class?????????? : RasterBrick
> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
> resolution : 1, 1?? (x, y)
> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> source???????? : memory
> names?????????? : layer.1, layer.2, layer.3
> min values :???????? 1.5,???????? 3.5,???????? 5.5
> max values :???????? 1.5,???????? 3.5,???????? 5.5??
>
>
> layer.1, layer.2, layer.3 (?)
>
> So what corrensponds to what?
>
>
> If I run:
>
> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>
> The result is:
>
> > res2
> class      : RasterBrick
> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
> resolution : 1, 1  (x, y)
> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> source     : memory
> names      : index_2, index_3, index_1
> min values :     1.5,     3.5,     5.5
> max values :     1.5,     3.5,     5.5
>
> There is no consistency with the names of the output and obscure
> correspondence with the indices in the case of clusterR
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
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 4:15 PM Leonidas Liakos via R-sig-Geo <
> [hidden email]> wrote:
>
>> I run the example with clusterR:
>>
>> no_cores <- parallel::detectCores() -1
>> raster::beginCluster(no_cores)
>> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>> list(indices=c(2,2,3,3,1,1),fun = mean))
>> raster::endCluster()
>>
>> And the result is:
>>
>>> res
>> class?????????? : RasterBrick
>> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>> resolution : 1, 1?? (x, y)
>> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> source???????? : memory
>> names?????????? : layer.1, layer.2, layer.3
>> min values :???????? 1.5,???????? 3.5,???????? 5.5
>> max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>
>>
>> layer.1, layer.2, layer.3 (?)
>>
>> So what corrensponds to what?
>>
>>
>> If I run:
>>
>> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>
>> The result is:
>>
>>> res2
>> class      : RasterBrick
>> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>> resolution : 1, 1  (x, y)
>> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> source     : memory
>> names      : index_2, index_3, index_1
>> min values :     1.5,     3.5,     5.5
>> max values :     1.5,     3.5,     5.5
>>
>> There is no consistency with the names of the output and obscure
>> correspondence with the indices in the case of clusterR
>>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>


--
Λιάκος Λεωνίδας, Γεωγράφος
https://www.geographer.gr
PGP fingerprint: 5237 83F8 E46C D91A 9FBB C7E7 F943 C9B6 8231 0937

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

Ben Tupper
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
<[hidden email]> 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 4:15 PM Leonidas Liakos via R-sig-Geo <
> > [hidden email]> wrote:
> >
> >> I run the example with clusterR:
> >>
> >> no_cores <- parallel::detectCores() -1
> >> raster::beginCluster(no_cores)
> >> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
> >> list(indices=c(2,2,3,3,1,1),fun = mean))
> >> raster::endCluster()
> >>
> >> And the result is:
> >>
> >>> res
> >> class?????????? : RasterBrick
> >> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
> >> resolution : 1, 1?? (x, y)
> >> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
> >> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >> source???????? : memory
> >> names?????????? : layer.1, layer.2, layer.3
> >> min values :???????? 1.5,???????? 3.5,???????? 5.5
> >> max values :???????? 1.5,???????? 3.5,???????? 5.5??
> >>
> >>
> >> layer.1, layer.2, layer.3 (?)
> >>
> >> So what corrensponds to what?
> >>
> >>
> >> If I run:
> >>
> >> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
> >>
> >> The result is:
> >>
> >>> res2
> >> class      : RasterBrick
> >> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
> >> resolution : 1, 1  (x, y)
> >> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> >> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >> source     : memory
> >> names      : index_2, index_3, index_1
> >> min values :     1.5,     3.5,     5.5
> >> max values :     1.5,     3.5,     5.5
> >>
> >> There is no consistency with the names of the output and obscure
> >> correspondence with the indices in the case of clusterR
> >>
> >>
> >>         [[alternative HTML version deleted]]
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
>
>
> --
> Λιάκος Λεωνίδας, Γεωγράφος
> https://www.geographer.gr
> PGP fingerprint: 5237 83F8 E46C D91A 9FBB C7E7 F943 C9B6 8231 0937
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



--
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
http://www.bigelow.org/
https://eco.bigelow.org

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
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://github.com/rspatial/raster/issues/82


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
> <[hidden email]> 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 4:15 PM Leonidas Liakos via R-sig-Geo <
>>> [hidden email]> wrote:
>>>
>>>> I run the example with clusterR:
>>>>
>>>> no_cores <- parallel::detectCores() -1
>>>> raster::beginCluster(no_cores)
>>>> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>>>> list(indices=c(2,2,3,3,1,1),fun = mean))
>>>> raster::endCluster()
>>>>
>>>> And the result is:
>>>>
>>>>> res
>>>> class?????????? : RasterBrick
>>>> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>>>> resolution : 1, 1?? (x, y)
>>>> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>>>> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>>>> source???????? : memory
>>>> names?????????? : layer.1, layer.2, layer.3
>>>> min values :???????? 1.5,???????? 3.5,???????? 5.5
>>>> max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>>>
>>>>
>>>> layer.1, layer.2, layer.3 (?)
>>>>
>>>> So what corrensponds to what?
>>>>
>>>>
>>>> If I run:
>>>>
>>>> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>>>
>>>> The result is:
>>>>
>>>>> res2
>>>> class      : RasterBrick
>>>> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>>>> resolution : 1, 1  (x, y)
>>>> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>>>> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>>>> source     : memory
>>>> names      : index_2, index_3, index_1
>>>> min values :     1.5,     3.5,     5.5
>>>> max values :     1.5,     3.5,     5.5
>>>>
>>>> There is no consistency with the names of the output and obscure
>>>> correspondence with the indices in the case of clusterR
>>>>
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>
>> --
>> Λιάκος Λεωνίδας, Γεωγράφος
>> https://www.geographer.gr
>> PGP fingerprint: 5237 83F8 E46C D91A 9FBB C7E7 F943 C9B6 8231 0937
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: 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
[hidden email]<https://remi.webmail.ec.europa.eu/owa/redir.aspx?C=O12RUARdbvGA3WF3zGoSV0j5xMoZlQcIEwiS4Y9G8jzXRqCCC1HUCA..&URL=mailto%3ajon.skoien%40jrc.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 <[hidden email]> on behalf of Leonidas Liakos via R-sig-Geo <[hidden email]>
Sent: 21 November 2019 08:52
To: Ben Tupper; [hidden email]
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
> <[hidden email]> 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 4:15 PM Leonidas Liakos via R-sig-Geo <
>>> [hidden email]> wrote:
>>>
>>>> I run the example with clusterR:
>>>>
>>>> no_cores <- parallel::detectCores() -1
>>>> raster::beginCluster(no_cores)
>>>> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>>>> list(indices=c(2,2,3,3,1,1),fun = mean))
>>>> raster::endCluster()
>>>>
>>>> And the result is:
>>>>
>>>>> res
>>>> class?????????? : RasterBrick
>>>> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>>>> resolution : 1, 1?? (x, y)
>>>> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>>>> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>>>> source???????? : memory
>>>> names?????????? : layer.1, layer.2, layer.3
>>>> min values :???????? 1.5,???????? 3.5,???????? 5.5
>>>> max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>>>
>>>>
>>>> layer.1, layer.2, layer.3 (?)
>>>>
>>>> So what corrensponds to what?
>>>>
>>>>
>>>> If I run:
>>>>
>>>> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>>>
>>>> The result is:
>>>>
>>>>> res2
>>>> class      : RasterBrick
>>>> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>>>> resolution : 1, 1  (x, y)
>>>> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>>>> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>>>> source     : memory
>>>> names      : index_2, index_3, index_1
>>>> min values :     1.5,     3.5,     5.5
>>>> max values :     1.5,     3.5,     5.5
>>>>
>>>> There is no consistency with the names of the output and obscure
>>>> correspondence with the indices in the case of clusterR
>>>>
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>>>>
>>
>> --
>> Λιάκος Λεωνίδας, Γεωγράφος
>> https://urldefense.com/v3/__https://www.geographer.gr__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835GzqxUtB7$
>> PGP fingerprint: 5237 83F8 E46C D91A 9FBB C7E7 F943 C9B6 8231 0937
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
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, 90000, 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=500000 +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, 90000, 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=500000 +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, [hidden email] 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
> [hidden email] 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 <[hidden email]> on behalf of
> Leonidas Liakos via R-sig-Geo <[hidden email]>
> *Sent:* 21 November 2019 08:52
> *To:* Ben Tupper; [hidden email]
> *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
> > <[hidden email]> 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 4:15 PM Leonidas Liakos via R-sig-Geo <
> >>> [hidden email]> wrote:
> >>>
> >>>> I run the example with clusterR:
> >>>>
> >>>> no_cores <- parallel::detectCores() -1
> >>>> raster::beginCluster(no_cores)
> >>>> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
> >>>> list(indices=c(2,2,3,3,1,1),fun = mean))
> >>>> raster::endCluster()
> >>>>
> >>>> And the result is:
> >>>>
> >>>>> res
> >>>> class?????????? : RasterBrick
> >>>> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
> >>>> resolution : 1, 1?? (x, y)
> >>>> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
> >>>> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
> +towgs84=0,0,0
> >>>> source???????? : memory
> >>>> names?????????? : layer.1, layer.2, layer.3
> >>>> min values :???????? 1.5,???????? 3.5,???????? 5.5
> >>>> max values :???????? 1.5,???????? 3.5,???????? 5.5??
> >>>>
> >>>>
> >>>> layer.1, layer.2, layer.3 (?)
> >>>>
> >>>> So what corrensponds to what?
> >>>>
> >>>>
> >>>> If I run:
> >>>>
> >>>> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
> >>>>
> >>>> The result is:
> >>>>
> >>>>> res2
> >>>> class      : RasterBrick
> >>>> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
> >>>> resolution : 1, 1  (x, y)
> >>>> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> >>>> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >>>> source     : memory
> >>>> names      : index_2, index_3, index_1
> >>>> min values :     1.5,     3.5,     5.5
> >>>> max values :     1.5,     3.5,     5.5
> >>>>
> >>>> There is no consistency with the names of the output and obscure
> >>>> correspondence with the indices in the case of clusterR
> >>>>
> >>>>
> >>>>         [[alternative HTML version deleted]]
> >>>>
> >>>> _______________________________________________
> >>>> R-sig-Geo mailing list
> >>>> [hidden email]
> >>>>
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
> >>>>
> >>
> >> --
> >> Λιάκος Λεωνίδας, Γεωγράφος
> >>
> https://urldefense.com/v3/__https://www.geographer.gr__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835GzqxUtB7$
>
> >> PGP fingerprint: 5237 83F8 E46C D91A 9FBB C7E7 F943 C9B6 8231 0937
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >>
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
> >
> >
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
In reply to this post by Jon.SKOIEN
Ok, I think it's a bug and I'd like your help to confirm this. When the
x parameter (a raster stack in this case) of stackApply is quite heavy,
then the result of the stackApply is written to the disk and for some
reason the name indexes are shuffled. Conversely, when the provided
raster stack  is lightweight, the result is stored in memory, and the
name indexes match the provided indices correctly. I am posting a gist
again where I create a heavy raster stack to replicate the problem.
Please ensure that the raster stack is heavy enough for your machine
that the result of the stackApply is written to the disk.
Could you please confirm it?
Here is the gist:
https://gist.github.com/kokkytos/76a283a770df10c05b4fb4ce97e2b40e



On 11/22/19 6:31 PM, [hidden email] 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
> [hidden email] 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 <[hidden email]> on behalf of
> Leonidas Liakos via R-sig-Geo <[hidden email]>
> *Sent:* 21 November 2019 08:52
> *To:* Ben Tupper; [hidden email]
> *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
> > <[hidden email]> 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 4:15 PM Leonidas Liakos via R-sig-Geo <
> >>> [hidden email]> wrote:
> >>>
> >>>> I run the example with clusterR:
> >>>>
> >>>> no_cores <- parallel::detectCores() -1
> >>>> raster::beginCluster(no_cores)
> >>>> ?????? res <- raster::clusterR(inp, raster::stackApply, args =
> >>>> list(indices=c(2,2,3,3,1,1),fun = mean))
> >>>> raster::endCluster()
> >>>>
> >>>> And the result is:
> >>>>
> >>>>> res
> >>>> class?????????? : RasterBrick
> >>>> dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
> >>>> resolution : 1, 1?? (x, y)
> >>>> extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
> >>>> crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
> +towgs84=0,0,0
> >>>> source???????? : memory
> >>>> names?????????? : layer.1, layer.2, layer.3
> >>>> min values :???????? 1.5,???????? 3.5,???????? 5.5
> >>>> max values :???????? 1.5,???????? 3.5,???????? 5.5??
> >>>>
> >>>>
> >>>> layer.1, layer.2, layer.3 (?)
> >>>>
> >>>> So what corrensponds to what?
> >>>>
> >>>>
> >>>> If I run:
> >>>>
> >>>> res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
> >>>>
> >>>> The result is:
> >>>>
> >>>>> res2
> >>>> class      : RasterBrick
> >>>> dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
> >>>> resolution : 1, 1  (x, y)
> >>>> extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> >>>> crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >>>> source     : memory
> >>>> names      : index_2, index_3, index_1
> >>>> min values :     1.5,     3.5,     5.5
> >>>> max values :     1.5,     3.5,     5.5
> >>>>
> >>>> There is no consistency with the names of the output and obscure
> >>>> correspondence with the indices in the case of clusterR
> >>>>
> >>>>
> >>>>         [[alternative HTML version deleted]]
> >>>>
> >>>> _______________________________________________
> >>>> R-sig-Geo mailing list
> >>>> [hidden email]
> >>>>
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
> >>>>
> >>
> >> --
> >> Λιάκος Λεωνίδας, Γεωγράφος
> >>
> https://urldefense.com/v3/__https://www.geographer.gr__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835GzqxUtB7$
>
> >> PGP fingerprint: 5237 83F8 E46C D91A 9FBB C7E7 F943 C9B6 8231 0937
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >>
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
> >
> >
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
In reply to this post by Frederico Faleiro
I added raster::zApply in my tests to validate the results. However, the
indices of the names of the results are different now. Recall that the
goal is to calculate from a raster stack time series the mean per day of
the week. And that problem I have is that stackApply, zApply and
calc/sapply return different indices in the result names. New code is
available here:
https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
I'm really curious about missing something.


On 11/20/19 3:30 AM, Frederico Faleiro wrote:

> 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 4:15 PM Leonidas Liakos via R-sig-Geo
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     I run the example with clusterR:
>
>     no_cores <- parallel::detectCores() -1
>     raster::beginCluster(no_cores)
>     ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>     list(indices=c(2,2,3,3,1,1),fun = mean))
>     raster::endCluster()
>
>     And the result is:
>
>     > res
>     class?????????? : RasterBrick
>     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>     resolution : 1, 1?? (x, y)
>     extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
>     +towgs84=0,0,0
>     source???????? : memory
>     names?????????? : layer.1, layer.2, layer.3
>     min values :???????? 1.5,???????? 3.5,???????? 5.5
>     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>
>
>     layer.1, layer.2, layer.3 (?)
>
>     So what corrensponds to what?
>
>
>     If I run:
>
>     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>
>     The result is:
>
>     > res2
>     class      : RasterBrick
>     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>     resolution : 1, 1  (x, y)
>     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>     source     : memory
>     names      : index_2, index_3, index_1
>     min values :     1.5,     3.5,     5.5
>     max values :     1.5,     3.5,     5.5
>
>     There is no consistency with the names of the output and obscure
>     correspondence with the indices in the case of clusterR
>
>
>             [[alternative HTML version deleted]]
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     [hidden email] <mailto:[hidden email]>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

Vijay Lulla
If you read the code/help for `stackApply` and `zApply` you'll see that the
results that you obtain make sense (at least they seem sensible/reasonable
to me).  IMO, if you want to control the ordering of your layers then just
use sapply, like how you've used for ver_mean.  IMO, this is the only
reliable (safe?), and quite a readable, way to accomplish what you're
trying to do.
Just my 2 cents.
-- Vijay.

On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via R-sig-Geo <
[hidden email]> wrote:

> I added raster::zApply in my tests to validate the results. However, the
> indices of the names of the results are different now. Recall that the
> goal is to calculate from a raster stack time series the mean per day of
> the week. And that problem I have is that stackApply, zApply and
> calc/sapply return different indices in the result names. New code is
> available here:
> https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
> I'm really curious about missing something.
>
>
> On 11/20/19 3:30 AM, Frederico Faleiro wrote:
> > 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 4:15 PM Leonidas Liakos via R-sig-Geo
> > <[hidden email] <mailto:[hidden email]>> wrote:
> >
> >     I run the example with clusterR:
> >
> >     no_cores <- parallel::detectCores() -1
> >     raster::beginCluster(no_cores)
> >     ?????? res <- raster::clusterR(inp, raster::stackApply, args =
> >     list(indices=c(2,2,3,3,1,1),fun = mean))
> >     raster::endCluster()
> >
> >     And the result is:
> >
> >     > res
> >     class?????????? : RasterBrick
> >     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
> >     resolution : 1, 1?? (x, y)
> >     extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
> >     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
> >     +towgs84=0,0,0
> >     source???????? : memory
> >     names?????????? : layer.1, layer.2, layer.3
> >     min values :???????? 1.5,???????? 3.5,???????? 5.5
> >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
> >
> >
> >     layer.1, layer.2, layer.3 (?)
> >
> >     So what corrensponds to what?
> >
> >
> >     If I run:
> >
> >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
> >
> >     The result is:
> >
> >     > res2
> >     class      : RasterBrick
> >     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
> >     resolution : 1, 1  (x, y)
> >     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> >     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >     source     : memory
> >     names      : index_2, index_3, index_1
> >     min values :     1.5,     3.5,     5.5
> >     max values :     1.5,     3.5,     5.5
> >
> >     There is no consistency with the names of the output and obscure
> >     correspondence with the indices in the case of clusterR
> >
> >
> >             [[alternative HTML version deleted]]
> >
> >     _______________________________________________
> >     R-sig-Geo mailing list
> >     [hidden email] <mailto:[hidden email]>
> >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
Why do they seem logical since they do not match?

Check for example index 1 (Sunday). The results are different for the
three processes

> stackapply_mean
class      : RasterBrick
dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
resolution : 500, 500  (x, y)
extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
crs        : NA
source     : /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
names      :  index_5,  index_6,  index_7,  index_1,  index_2, 
index_3,  index_4
min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
442.7436
max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
556.2471


> ver_mean
class      : RasterStack
dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
resolution : 500, 500  (x, y)
extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
crs        : NA
names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5, 
layer.6,  layer.7
min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
429.6900
max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
561.7972


> z
class      : RasterBrick
dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
resolution : 500, 500  (x, y)
extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
crs        : NA
source     : /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
names      :       X1,       X2,       X3,       X4,       X5,      
X6,       X7
min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
442.7436
max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
556.2471
           : 1, 2, 3, 4, 5, 6, 7


On 11/26/19 7:03 PM, Vijay Lulla wrote:

> If you read the code/help for `stackApply` and `zApply` you'll see
> that the results that you obtain make sense (at least they seem
> sensible/reasonable to me).  IMO, if you want to control the ordering
> of your layers then just use sapply, like how you've used for
> ver_mean.  IMO, this is the only reliable (safe?), and quite a
> readable, way to accomplish what you're trying to do.
> Just my 2 cents.
> -- Vijay.
>
> On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via R-sig-Geo
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     I added raster::zApply in my tests to validate the results.
>     However, the
>     indices of the names of the results are different now. Recall that the
>     goal is to calculate from a raster stack time series the mean per
>     day of
>     the week. And that problem I have is that stackApply, zApply and
>     calc/sapply return different indices in the result names. New code is
>     available here:
>     https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>     I'm really curious about missing something.
>
>
>     On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>     > 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 4:15 PM Leonidas Liakos via R-sig-Geo
>     > <[hidden email] <mailto:[hidden email]>
>     <mailto:[hidden email] <mailto:[hidden email]>>>
>     wrote:
>     >
>     >     I run the example with clusterR:
>     >
>     >     no_cores <- parallel::detectCores() -1
>     >     raster::beginCluster(no_cores)
>     >     ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>     >     list(indices=c(2,2,3,3,1,1),fun = mean))
>     >     raster::endCluster()
>     >
>     >     And the result is:
>     >
>     >     > res
>     >     class?????????? : RasterBrick
>     >     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>     >     resolution : 1, 1?? (x, y)
>     >     extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>     >     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
>     >     +towgs84=0,0,0
>     >     source???????? : memory
>     >     names?????????? : layer.1, layer.2, layer.3
>     >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>     >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>     >
>     >
>     >     layer.1, layer.2, layer.3 (?)
>     >
>     >     So what corrensponds to what?
>     >
>     >
>     >     If I run:
>     >
>     >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>     >
>     >     The result is:
>     >
>     >     > res2
>     >     class      : RasterBrick
>     >     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>     >     resolution : 1, 1  (x, y)
>     >     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>     >     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84
>     +towgs84=0,0,0
>     >     source     : memory
>     >     names      : index_2, index_3, index_1
>     >     min values :     1.5,     3.5,     5.5
>     >     max values :     1.5,     3.5,     5.5
>     >
>     >     There is no consistency with the names of the output and obscure
>     >     correspondence with the indices in the case of clusterR
>     >
>     >
>     >             [[alternative HTML version deleted]]
>     >
>     >     _______________________________________________
>     >     R-sig-Geo mailing list
>     >     [hidden email] <mailto:[hidden email]>
>     <mailto:[hidden email] <mailto:[hidden email]>>
>     >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     >
>
>             [[alternative HTML version deleted]]
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     [hidden email] <mailto:[hidden email]>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

Vijay Lulla
I'm sorry for the miscommunication.  What I meant to say is that the output
from stackApply and zApply are the same (because zApply uses stackApply
internally) except the names.  The behavior of stackApply makes sense
because AFAIUI R doesn't automatically reorder vectors/indices that it
receives.  Your observation about inconsistent result with ver_mean is very
valid though!  And, that's what I meant with my comment that using sapply
with the explicit ordering that you want is the best way to control what
raster package will output.  In R the input order should be maintained
(this is the prime difference between SQL and R) but packages/tools do not
always adhere to this...so it's never clear how the output will be
ordered.  Sorry for the confusion.


On Tue, Nov 26, 2019 at 12:22 PM Leonidas Liakos <[hidden email]>
wrote:

> Why do they seem logical since they do not match?
>
> Check for example index 1 (Sunday). The results are different for the
> three processes
>
> > stackapply_mean
> class      : RasterBrick
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> source     : /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
> names      :  index_5,  index_6,  index_7,  index_1,  index_2,  index_3,
> index_4
> min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
> 442.7436
> max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
> 556.2471
>
>
> > ver_mean
> class      : RasterStack
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,
> layer.7
> min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
> 429.6900
> max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
> 561.7972
>
>
> > z
> class      : RasterBrick
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> source     : /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
> names      :       X1,       X2,       X3,       X4,       X5,
> X6,       X7
> min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
> 442.7436
> max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
> 556.2471
>            : 1, 2, 3, 4, 5, 6, 7
>
>
> On 11/26/19 7:03 PM, Vijay Lulla wrote:
>
> If you read the code/help for `stackApply` and `zApply` you'll see that
> the results that you obtain make sense (at least they seem
> sensible/reasonable to me).  IMO, if you want to control the ordering of
> your layers then just use sapply, like how you've used for ver_mean.  IMO,
> this is the only reliable (safe?), and quite a readable, way to accomplish
> what you're trying to do.
> Just my 2 cents.
> -- Vijay.
>
> On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via R-sig-Geo <
> [hidden email]> wrote:
>
>> I added raster::zApply in my tests to validate the results. However, the
>> indices of the names of the results are different now. Recall that the
>> goal is to calculate from a raster stack time series the mean per day of
>> the week. And that problem I have is that stackApply, zApply and
>> calc/sapply return different indices in the result names. New code is
>> available here:
>> https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>> I'm really curious about missing something.
>>
>>
>> On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>> > 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 4:15 PM Leonidas Liakos via R-sig-Geo
>> > <[hidden email] <mailto:[hidden email]>> wrote:
>> >
>> >     I run the example with clusterR:
>> >
>> >     no_cores <- parallel::detectCores() -1
>> >     raster::beginCluster(no_cores)
>> >     ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>> >     list(indices=c(2,2,3,3,1,1),fun = mean))
>> >     raster::endCluster()
>> >
>> >     And the result is:
>> >
>> >     > res
>> >     class?????????? : RasterBrick
>> >     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>> >     resolution : 1, 1?? (x, y)
>> >     extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>> >     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
>> >     +towgs84=0,0,0
>> >     source???????? : memory
>> >     names?????????? : layer.1, layer.2, layer.3
>> >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>> >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>> >
>> >
>> >     layer.1, layer.2, layer.3 (?)
>> >
>> >     So what corrensponds to what?
>> >
>> >
>> >     If I run:
>> >
>> >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>> >
>> >     The result is:
>> >
>> >     > res2
>> >     class      : RasterBrick
>> >     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>> >     resolution : 1, 1  (x, y)
>> >     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>> >     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> >     source     : memory
>> >     names      : index_2, index_3, index_1
>> >     min values :     1.5,     3.5,     5.5
>> >     max values :     1.5,     3.5,     5.5
>> >
>> >     There is no consistency with the names of the output and obscure
>> >     correspondence with the indices in the case of clusterR
>> >
>> >
>> >             [[alternative HTML version deleted]]
>> >
>> >     _______________________________________________
>> >     R-sig-Geo mailing list
>> >     [hidden email] <mailto:[hidden email]>
>> >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
Thank you!
The problem is not with the resulting values ​​but with the index
mapping. Values ​​are correct in all three cases.

As I wrote in a previous post in the thread
(https://stat.ethz.ch/pipermail/r-sig-geo/2019-November/027821.html) ,
stackApply behaves inconsistently depending on whether the exported
stack will remain in memory or it will be stored, due to its large size,
on the hard disk.

In the first case the indices are identical to my function (ver_mean)
and the lubridate::wday indexing system (and are correct) while in the
second they are shuffled.

So, while Sunday has index 1 and while in the first case (when the
result is in memory) stackApply returns the correct index, in the second
case (when the result is stored on the hard disk) it returns index_4! So
how can one be sure if index_1 corresponds to Sunday or another day
using stackApply since it sometimes enumerates it with index_1 and
sometimes index_4?


Try to run this example (when the resulting stack remains in memory) to
see that the indexes are identical (stackApply = ver_median =
lubridate::wday)
https://gist.github.com/kokkytos/5d554b5a725bb48d2189e2d1fa0e2206

Thank you again

On 11/26/19 9:00 PM, Vijay Lulla wrote:

> I'm sorry for the miscommunication.  What I meant to say is that the
> output from stackApply and zApply are the same (because zApply uses
> stackApply internally) except the names.  The behavior of stackApply
> makes sense because AFAIUI R doesn't automatically reorder
> vectors/indices that it receives.  Your observation about inconsistent
> result with ver_mean is very valid though!  And, that's what I meant
> with my comment that using sapply with the explicit ordering that you
> want is the best way to control what raster package will output.  In R
> the input order should be maintained (this is the prime difference
> between SQL and R) but packages/tools do not always adhere to
> this...so it's never clear how the output will be ordered.  Sorry for
> the confusion.
>
>
> On Tue, Nov 26, 2019 at 12:22 PM Leonidas Liakos
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     Why do they seem logical since they do not match?
>
>     Check for example index 1 (Sunday). The results are different for
>     the three processes
>
>     > stackapply_mean
>     class      : RasterBrick
>     dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>     resolution : 500, 500  (x, y)
>     extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>     crs        : NA
>     source     :
>     /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
>     names      :  index_5,  index_6,  index_7,  index_1,  index_2, 
>     index_3,  index_4
>     min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
>     429.6900, 442.7436
>     max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
>     561.7972, 556.2471
>
>
>     > ver_mean
>     class      : RasterStack
>     dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>     resolution : 500, 500  (x, y)
>     extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>     crs        : NA
>     names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5, 
>     layer.6,  layer.7
>     min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946,
>     440.2028, 429.6900
>     max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671,
>     560.1375, 561.7972
>
>
>     > z
>     class      : RasterBrick
>     dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>     resolution : 500, 500  (x, y)
>     extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>     crs        : NA
>     source     :
>     /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
>     names      :       X1,       X2,       X3,       X4,      
>     X5,       X6,       X7
>     min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
>     429.6900, 442.7436
>     max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
>     561.7972, 556.2471
>                : 1, 2, 3, 4, 5, 6, 7
>
>
>     On 11/26/19 7:03 PM, Vijay Lulla wrote:
>>     If you read the code/help for `stackApply` and `zApply` you'll
>>     see that the results that you obtain make sense (at least they
>>     seem sensible/reasonable to me).  IMO, if you want to control the
>>     ordering of your layers then just use sapply, like how you've
>>     used for ver_mean.  IMO, this is the only reliable (safe?), and
>>     quite a readable, way to accomplish what you're trying to do.
>>     Just my 2 cents.
>>     -- Vijay.
>>
>>     On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via R-sig-Geo
>>     <[hidden email] <mailto:[hidden email]>> wrote:
>>
>>         I added raster::zApply in my tests to validate the results.
>>         However, the
>>         indices of the names of the results are different now. Recall
>>         that the
>>         goal is to calculate from a raster stack time series the mean
>>         per day of
>>         the week. And that problem I have is that stackApply, zApply and
>>         calc/sapply return different indices in the result names. New
>>         code is
>>         available here:
>>         https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>>         I'm really curious about missing something.
>>
>>
>>         On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>>         > 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 4:15 PM Leonidas Liakos via R-sig-Geo
>>         > <[hidden email] <mailto:[hidden email]>
>>         <mailto:[hidden email]
>>         <mailto:[hidden email]>>> wrote:
>>         >
>>         >     I run the example with clusterR:
>>         >
>>         >     no_cores <- parallel::detectCores() -1
>>         >     raster::beginCluster(no_cores)
>>         >     ?????? res <- raster::clusterR(inp, raster::stackApply,
>>         args =
>>         >     list(indices=c(2,2,3,3,1,1),fun = mean))
>>         >     raster::endCluster()
>>         >
>>         >     And the result is:
>>         >
>>         >     > res
>>         >     class?????????? : RasterBrick
>>         >     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell,
>>         nlayers)
>>         >     resolution : 1, 1?? (x, y)
>>         >     extent???????? : -180, 180, -90, 90?? (xmin, xmax,
>>         ymin, ymax)
>>         >     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
>>         >     +towgs84=0,0,0
>>         >     source???????? : memory
>>         >     names?????????? : layer.1, layer.2, layer.3
>>         >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>>         >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>         >
>>         >
>>         >     layer.1, layer.2, layer.3 (?)
>>         >
>>         >     So what corrensponds to what?
>>         >
>>         >
>>         >     If I run:
>>         >
>>         >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>         >
>>         >     The result is:
>>         >
>>         >     > res2
>>         >     class      : RasterBrick
>>         >     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell,
>>         nlayers)
>>         >     resolution : 1, 1  (x, y)
>>         >     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>>         >     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84
>>         +towgs84=0,0,0
>>         >     source     : memory
>>         >     names      : index_2, index_3, index_1
>>         >     min values :     1.5,     3.5,     5.5
>>         >     max values :     1.5,     3.5,     5.5
>>         >
>>         >     There is no consistency with the names of the output
>>         and obscure
>>         >     correspondence with the indices in the case of clusterR
>>         >
>>         >
>>         >             [[alternative HTML version deleted]]
>>         >
>>         >     _______________________________________________
>>         >     R-sig-Geo mailing list
>>         >     [hidden email]
>>         <mailto:[hidden email]>
>>         <mailto:[hidden email] <mailto:[hidden email]>>
>>         >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>         >
>>
>>                 [[alternative HTML version deleted]]
>>
>>         _______________________________________________
>>         R-sig-Geo mailing list
>>         [hidden email] <mailto:[hidden email]>
>>         https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

Vijay Lulla
Hmm...it appears that stackApply is using different conditions which might
be causing this problem. Below is the snippet of the code which I think
might be the problem.

## For canProcessInMemory
if (rowcalc) {
  v <- lapply(uin, function(i) fun(x[, ind == i, drop = FALSE], na.rm =
na.rm))
}
else {
  v <- lapply(uin, function(i, ...) apply(x[, ind == i, drop = FALSE], 1,
fun, na.rm = na.rm))
}


## If canProcessInMemory is not TRUE
if (rowcalc) {
  v <- lapply(uin, function(i) fun(a[, ind == uin[i], drop = FALSE], na.rm
= na.rm))
}
else {
  v <- lapply(uin, function(i, ...) apply(a[, ind == uin[i], drop = FALSE],
1, fun, na.rm = na.rm))
}

I think they should both be same but it appears that they aren't and that's
what you've discovered.  Maybe you can try fix(stackApply) to see if it
really is the problem and can tell us what you find.  Anyways, good
catch...and...sorry for wasting your time.
Cordially,
Vijay.

On Tue, Nov 26, 2019 at 2:53 PM Leonidas Liakos <[hidden email]>
wrote:

> Thank you!
> The problem is not with the resulting values but with the index mapping.
> Values are correct in all three cases.
>
> As I wrote in a previous post in the thread (
> https://stat.ethz.ch/pipermail/r-sig-geo/2019-November/027821.html) ,
> stackApply behaves inconsistently depending on whether the exported stack
> will remain in memory or it will be stored, due to its large size, on the
> hard disk.
>
> In the first case the indices are identical to my function (ver_mean) and
> the lubridate::wday indexing system (and are correct) while in the second
> they are shuffled.
>
> So, while Sunday has index 1 and while in the first case (when the result
> is in memory) stackApply returns the correct index, in the second case
> (when the result is stored on the hard disk) it returns index_4! So how can
> one be sure if index_1 corresponds to Sunday or another day using
> stackApply since it sometimes enumerates it with index_1 and sometimes
> index_4?
>
>
> Try to run this example (when the resulting stack remains in memory) to
> see that the indexes are identical (stackApply = ver_median =
> lubridate::wday)
> https://gist.github.com/kokkytos/5d554b5a725bb48d2189e2d1fa0e2206
>
> Thank you again
> On 11/26/19 9:00 PM, Vijay Lulla wrote:
>
> I'm sorry for the miscommunication.  What I meant to say is that the
> output from stackApply and zApply are the same (because zApply uses
> stackApply internally) except the names.  The behavior of stackApply makes
> sense because AFAIUI R doesn't automatically reorder vectors/indices that
> it receives.  Your observation about inconsistent result with ver_mean is
> very valid though!  And, that's what I meant with my comment that using
> sapply with the explicit ordering that you want is the best way to control
> what raster package will output.  In R the input order should be maintained
> (this is the prime difference between SQL and R) but packages/tools do not
> always adhere to this...so it's never clear how the output will be
> ordered.  Sorry for the confusion.
>
>
> On Tue, Nov 26, 2019 at 12:22 PM Leonidas Liakos <[hidden email]>
> wrote:
>
>> Why do they seem logical since they do not match?
>>
>> Check for example index 1 (Sunday). The results are different for the
>> three processes
>>
>> > stackapply_mean
>> class      : RasterBrick
>> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>> resolution : 500, 500  (x, y)
>> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>> crs        : NA
>> source     :
>> /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
>> names      :  index_5,  index_6,  index_7,  index_1,  index_2,  index_3,
>> index_4
>> min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
>> 442.7436
>> max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
>> 556.2471
>>
>>
>> > ver_mean
>> class      : RasterStack
>> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>> resolution : 500, 500  (x, y)
>> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>> crs        : NA
>> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,
>> layer.7
>> min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
>> 429.6900
>> max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
>> 561.7972
>>
>>
>> > z
>> class      : RasterBrick
>> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>> resolution : 500, 500  (x, y)
>> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>> crs        : NA
>> source     :
>> /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
>> names      :       X1,       X2,       X3,       X4,       X5,
>> X6,       X7
>> min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
>> 442.7436
>> max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
>> 556.2471
>>            : 1, 2, 3, 4, 5, 6, 7
>>
>>
>> On 11/26/19 7:03 PM, Vijay Lulla wrote:
>>
>> If you read the code/help for `stackApply` and `zApply` you'll see that
>> the results that you obtain make sense (at least they seem
>> sensible/reasonable to me).  IMO, if you want to control the ordering of
>> your layers then just use sapply, like how you've used for ver_mean.  IMO,
>> this is the only reliable (safe?), and quite a readable, way to accomplish
>> what you're trying to do.
>> Just my 2 cents.
>> -- Vijay.
>>
>> On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via R-sig-Geo <
>> [hidden email]> wrote:
>>
>>> I added raster::zApply in my tests to validate the results. However, the
>>> indices of the names of the results are different now. Recall that the
>>> goal is to calculate from a raster stack time series the mean per day of
>>> the week. And that problem I have is that stackApply, zApply and
>>> calc/sapply return different indices in the result names. New code is
>>> available here:
>>> https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>>> I'm really curious about missing something.
>>>
>>>
>>> On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>>> > 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 4:15 PM Leonidas Liakos via R-sig-Geo
>>> > <[hidden email] <mailto:[hidden email]>> wrote:
>>> >
>>> >     I run the example with clusterR:
>>> >
>>> >     no_cores <- parallel::detectCores() -1
>>> >     raster::beginCluster(no_cores)
>>> >     ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>>> >     list(indices=c(2,2,3,3,1,1),fun = mean))
>>> >     raster::endCluster()
>>> >
>>> >     And the result is:
>>> >
>>> >     > res
>>> >     class?????????? : RasterBrick
>>> >     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>>> >     resolution : 1, 1?? (x, y)
>>> >     extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>>> >     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
>>> >     +towgs84=0,0,0
>>> >     source???????? : memory
>>> >     names?????????? : layer.1, layer.2, layer.3
>>> >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>>> >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>> >
>>> >
>>> >     layer.1, layer.2, layer.3 (?)
>>> >
>>> >     So what corrensponds to what?
>>> >
>>> >
>>> >     If I run:
>>> >
>>> >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>> >
>>> >     The result is:
>>> >
>>> >     > res2
>>> >     class      : RasterBrick
>>> >     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>>> >     resolution : 1, 1  (x, y)
>>> >     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>>> >     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>>> >     source     : memory
>>> >     names      : index_2, index_3, index_1
>>> >     min values :     1.5,     3.5,     5.5
>>> >     max values :     1.5,     3.5,     5.5
>>> >
>>> >     There is no consistency with the names of the output and obscure
>>> >     correspondence with the indices in the case of clusterR
>>> >
>>> >
>>> >             [[alternative HTML version deleted]]
>>> >
>>> >     _______________________________________________
>>> >     R-sig-Geo mailing list
>>> >     [hidden email] <mailto:[hidden email]>
>>> >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> >
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

R-sig-geo mailing list
Thank you for your help!

I tried to fix stackApply according to your instructions.

Now the indices of names are the same and consistent with indices
enumeration (gist for validation and tests:
https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170).

I've attached a patch file here:

https://gist.github.com/kokkytos/ca2c319134677b19900579665267a7a7

> stackapply_mean
class      : RasterBrick
dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
resolution : 500, 500  (x, y)
extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
crs        : NA
source     : /tmp/Rtmp9W8UNc/raster/r_tmp_2019-11-27_191205_2929_20324.grd
names      :  index_5,  index_6,  index_7,  index_1,  index_2, 
index_3,  index_4
min values : 444.6946, 440.2028, 429.6900, 442.7436, 440.0467, 444.9182,
437.1589
max values : 565.8671, 560.1375, 561.7972, 556.2471, 563.8341, 561.7687,
560.4509

> ver_mean
class      : RasterStack
dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
resolution : 500, 500  (x, y)
extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
crs        : NA
names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5, 
layer.6,  layer.7
min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
429.6900
max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
561.7972





On 11/26/19 10:58 PM, Vijay Lulla wrote:

> Hmm...it appears that stackApply is using different conditions which
> might be causing this problem. Below is the snippet of the code which
> I think might be the problem.
>
> ## For canProcessInMemory
> if (rowcalc) {
>   v <- lapply(uin, function(i) fun(x[, ind == i, drop = FALSE], na.rm
> = na.rm))
> }
> else {
>   v <- lapply(uin, function(i, ...) apply(x[, ind == i, drop = FALSE],
> 1, fun, na.rm = na.rm))
> }
>
>
> ## If canProcessInMemory is not TRUE
> if (rowcalc) {
>   v <- lapply(uin, function(i) fun(a[, ind == uin[i], drop = FALSE],
> na.rm = na.rm))
> }
> else {
>   v <- lapply(uin, function(i, ...) apply(a[, ind == uin[i], drop =
> FALSE], 1, fun, na.rm = na.rm))
> }
>
> I think they should both be same but it appears that they aren't and
> that's what you've discovered.  Maybe you can try fix(stackApply) to
> see if it really is the problem and can tell us what you find. 
> Anyways, good catch...and...sorry for wasting your time.
> Cordially,
> Vijay.
>
> On Tue, Nov 26, 2019 at 2:53 PM Leonidas Liakos
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     Thank you!
>     The problem is not with the resulting values but with the index
>     mapping. Values are correct in all three cases.
>
>     As I wrote in a previous post in the thread
>     (https://stat.ethz.ch/pipermail/r-sig-geo/2019-November/027821.html)
>     , stackApply behaves inconsistently depending on whether the
>     exported stack will remain in memory or it will be stored, due to
>     its large size, on the hard disk.
>
>     In the first case the indices are identical to my function
>     (ver_mean) and the lubridate::wday indexing system (and are
>     correct) while in the second they are shuffled.
>
>     So, while Sunday has index 1 and while in the first case (when the
>     result is in memory) stackApply returns the correct index, in the
>     second case (when the result is stored on the hard disk) it
>     returns index_4! So how can one be sure if index_1 corresponds to
>     Sunday or another day using stackApply since it sometimes
>     enumerates it with index_1 and sometimes index_4?
>
>
>     Try to run this example (when the resulting stack remains in
>     memory) to see that the indexes are identical (stackApply =
>     ver_median = lubridate::wday)
>     https://gist.github.com/kokkytos/5d554b5a725bb48d2189e2d1fa0e2206
>
>     Thank you again
>
>     On 11/26/19 9:00 PM, Vijay Lulla wrote:
>>     I'm sorry for the miscommunication.  What I meant to say is that
>>     the output from stackApply and zApply are the same (because
>>     zApply uses stackApply internally) except the names.  The
>>     behavior of stackApply makes sense because AFAIUI R doesn't
>>     automatically reorder vectors/indices that it receives.  Your
>>     observation about inconsistent result with ver_mean is very valid
>>     though!  And, that's what I meant with my comment that using
>>     sapply with the explicit ordering that you want is the best way
>>     to control what raster package will output.  In R the input order
>>     should be maintained (this is the prime difference between SQL
>>     and R) but packages/tools do not always adhere to this...so it's
>>     never clear how the output will be ordered.  Sorry for the confusion.
>>
>>
>>     On Tue, Nov 26, 2019 at 12:22 PM Leonidas Liakos
>>     <[hidden email] <mailto:[hidden email]>> wrote:
>>
>>         Why do they seem logical since they do not match?
>>
>>         Check for example index 1 (Sunday). The results are different
>>         for the three processes
>>
>>         > stackapply_mean
>>         class      : RasterBrick
>>         dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>         resolution : 500, 500  (x, y)
>>         extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>         crs        : NA
>>         source     :
>>         /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
>>         names      :  index_5,  index_6,  index_7,  index_1, 
>>         index_2,  index_3,  index_4
>>         min values : 440.0467, 444.9182, 437.1589, 444.6946,
>>         440.2028, 429.6900, 442.7436
>>         max values : 563.8341, 561.7687, 560.4509, 565.8671,
>>         560.1375, 561.7972, 556.2471
>>
>>
>>         > ver_mean
>>         class      : RasterStack
>>         dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>         resolution : 500, 500  (x, y)
>>         extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>         crs        : NA
>>         names      :  layer.1,  layer.2,  layer.3,  layer.4, 
>>         layer.5,  layer.6,  layer.7
>>         min values : 442.7436, 440.0467, 444.9182, 437.1589,
>>         444.6946, 440.2028, 429.6900
>>         max values : 556.2471, 563.8341, 561.7687, 560.4509,
>>         565.8671, 560.1375, 561.7972
>>
>>
>>         > z
>>         class      : RasterBrick
>>         dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>         resolution : 500, 500  (x, y)
>>         extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>         crs        : NA
>>         source     :
>>         /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
>>         names      :       X1,       X2,       X3,       X4,      
>>         X5,       X6,       X7
>>         min values : 440.0467, 444.9182, 437.1589, 444.6946,
>>         440.2028, 429.6900, 442.7436
>>         max values : 563.8341, 561.7687, 560.4509, 565.8671,
>>         560.1375, 561.7972, 556.2471
>>                    : 1, 2, 3, 4, 5, 6, 7
>>
>>
>>         On 11/26/19 7:03 PM, Vijay Lulla wrote:
>>>         If you read the code/help for `stackApply` and `zApply`
>>>         you'll see that the results that you obtain make sense (at
>>>         least they seem sensible/reasonable to me).  IMO, if you
>>>         want to control the ordering of your layers then just use
>>>         sapply, like how you've used for ver_mean.  IMO, this is the
>>>         only reliable (safe?), and quite a readable, way to
>>>         accomplish what you're trying to do.
>>>         Just my 2 cents.
>>>         -- Vijay.
>>>
>>>         On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via
>>>         R-sig-Geo <[hidden email]
>>>         <mailto:[hidden email]>> wrote:
>>>
>>>             I added raster::zApply in my tests to validate the
>>>             results. However, the
>>>             indices of the names of the results are different now.
>>>             Recall that the
>>>             goal is to calculate from a raster stack time series the
>>>             mean per day of
>>>             the week. And that problem I have is that stackApply,
>>>             zApply and
>>>             calc/sapply return different indices in the result
>>>             names. New code is
>>>             available here:
>>>             https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>>>             I'm really curious about missing something.
>>>
>>>
>>>             On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>>>             > 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 4:15 PM Leonidas Liakos via
>>>             R-sig-Geo
>>>             > <[hidden email]
>>>             <mailto:[hidden email]>
>>>             <mailto:[hidden email]
>>>             <mailto:[hidden email]>>> wrote:
>>>             >
>>>             >     I run the example with clusterR:
>>>             >
>>>             >     no_cores <- parallel::detectCores() -1
>>>             >     raster::beginCluster(no_cores)
>>>             >     ?????? res <- raster::clusterR(inp,
>>>             raster::stackApply, args =
>>>             >     list(indices=c(2,2,3,3,1,1),fun = mean))
>>>             >     raster::endCluster()
>>>             >
>>>             >     And the result is:
>>>             >
>>>             >     > res
>>>             >     class?????????? : RasterBrick
>>>             >     dimensions : 180, 360, 64800, 3?? (nrow, ncol,
>>>             ncell, nlayers)
>>>             >     resolution : 1, 1?? (x, y)
>>>             >     extent???????? : -180, 180, -90, 90?? (xmin, xmax,
>>>             ymin, ymax)
>>>             >     crs?????????????? : +proj=longlat +datum=WGS84
>>>             +ellps=WGS84
>>>             >     +towgs84=0,0,0
>>>             >     source???????? : memory
>>>             >     names?????????? : layer.1, layer.2, layer.3
>>>             >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>>>             >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>>             >
>>>             >
>>>             >     layer.1, layer.2, layer.3 (?)
>>>             >
>>>             >     So what corrensponds to what?
>>>             >
>>>             >
>>>             >     If I run:
>>>             >
>>>             >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>>             >
>>>             >     The result is:
>>>             >
>>>             >     > res2
>>>             >     class      : RasterBrick
>>>             >     dimensions : 180, 360, 64800, 3  (nrow, ncol,
>>>             ncell, nlayers)
>>>             >     resolution : 1, 1  (x, y)
>>>             >     extent     : -180, 180, -90, 90  (xmin, xmax,
>>>             ymin, ymax)
>>>             >     crs        : +proj=longlat +datum=WGS84
>>>             +ellps=WGS84 +towgs84=0,0,0
>>>             >     source     : memory
>>>             >     names      : index_2, index_3, index_1
>>>             >     min values :     1.5,     3.5,     5.5
>>>             >     max values :     1.5,     3.5,     5.5
>>>             >
>>>             >     There is no consistency with the names of the
>>>             output and obscure
>>>             >     correspondence with the indices in the case of
>>>             clusterR
>>>             >
>>>             >
>>>             >             [[alternative HTML version deleted]]
>>>             >
>>>             >     _______________________________________________
>>>             >     R-sig-Geo mailing list
>>>             >     [hidden email]
>>>             <mailto:[hidden email]>
>>>             <mailto:[hidden email]
>>>             <mailto:[hidden email]>>
>>>             >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>             >
>>>
>>>                     [[alternative HTML version deleted]]
>>>
>>>             _______________________________________________
>>>             R-sig-Geo mailing list
>>>             [hidden email] <mailto:[hidden email]>
>>>             https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>>
>>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

Vijay Lulla
Very nice!

And, just fyi: if you wish to get stackapply_mean to have layers in same
order as ver_mean then you can use  stackapply_mean <-
stack(stackapply_mean, layers=order(names(stackapply_mean)))

On Wed, Nov 27, 2019 at 12:43 PM Leonidas Liakos <[hidden email]>
wrote:

> Thank you for your help!
>
> I tried to fix stackApply according to your instructions.
>
> Now the indices of names are the same and consistent with indices
> enumeration (gist for validation and tests:
> https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170).
>
> I've attached a patch file here:
>
> https://gist.github.com/kokkytos/ca2c319134677b19900579665267a7a7
> > stackapply_mean
> class      : RasterBrick
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> source     : /tmp/Rtmp9W8UNc/raster/r_tmp_2019-11-27_191205_2929_20324.grd
> names      :  index_5,  index_6,  index_7,  index_1,  index_2,  index_3,
> index_4
> min values : 444.6946, 440.2028, 429.6900, 442.7436, 440.0467, 444.9182,
> 437.1589
> max values : 565.8671, 560.1375, 561.7972, 556.2471, 563.8341, 561.7687,
> 560.4509
>
> > ver_mean
> class      : RasterStack
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,
> layer.7
> min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
> 429.6900
> max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
> 561.7972
>
>
>
>
>
> On 11/26/19 10:58 PM, Vijay Lulla wrote:
>
> Hmm...it appears that stackApply is using different conditions which might
> be causing this problem. Below is the snippet of the code which I think
> might be the problem.
>
> ## For canProcessInMemory
> if (rowcalc) {
>   v <- lapply(uin, function(i) fun(x[, ind == i, drop = FALSE], na.rm =
> na.rm))
> }
> else {
>   v <- lapply(uin, function(i, ...) apply(x[, ind == i, drop = FALSE], 1,
> fun, na.rm = na.rm))
> }
>
>
> ## If canProcessInMemory is not TRUE
> if (rowcalc) {
>   v <- lapply(uin, function(i) fun(a[, ind == uin[i], drop = FALSE], na.rm
> = na.rm))
> }
> else {
>   v <- lapply(uin, function(i, ...) apply(a[, ind == uin[i], drop =
> FALSE], 1, fun, na.rm = na.rm))
> }
>
> I think they should both be same but it appears that they aren't and
> that's what you've discovered.  Maybe you can try fix(stackApply) to see if
> it really is the problem and can tell us what you find.  Anyways, good
> catch...and...sorry for wasting your time.
> Cordially,
> Vijay.
>
> On Tue, Nov 26, 2019 at 2:53 PM Leonidas Liakos <[hidden email]>
> wrote:
>
>> Thank you!
>> The problem is not with the resulting values but with the index mapping.
>> Values are correct in all three cases.
>>
>> As I wrote in a previous post in the thread (
>> https://stat.ethz.ch/pipermail/r-sig-geo/2019-November/027821.html) ,
>> stackApply behaves inconsistently depending on whether the exported stack
>> will remain in memory or it will be stored, due to its large size, on the
>> hard disk.
>>
>> In the first case the indices are identical to my function (ver_mean) and
>> the lubridate::wday indexing system (and are correct) while in the second
>> they are shuffled.
>>
>> So, while Sunday has index 1 and while in the first case (when the result
>> is in memory) stackApply returns the correct index, in the second case
>> (when the result is stored on the hard disk) it returns index_4! So how can
>> one be sure if index_1 corresponds to Sunday or another day using
>> stackApply since it sometimes enumerates it with index_1 and sometimes
>> index_4?
>>
>>
>> Try to run this example (when the resulting stack remains in memory) to
>> see that the indexes are identical (stackApply = ver_median =
>> lubridate::wday)
>> https://gist.github.com/kokkytos/5d554b5a725bb48d2189e2d1fa0e2206
>>
>> Thank you again
>> On 11/26/19 9:00 PM, Vijay Lulla wrote:
>>
>> I'm sorry for the miscommunication.  What I meant to say is that the
>> output from stackApply and zApply are the same (because zApply uses
>> stackApply internally) except the names.  The behavior of stackApply makes
>> sense because AFAIUI R doesn't automatically reorder vectors/indices that
>> it receives.  Your observation about inconsistent result with ver_mean is
>> very valid though!  And, that's what I meant with my comment that using
>> sapply with the explicit ordering that you want is the best way to control
>> what raster package will output.  In R the input order should be maintained
>> (this is the prime difference between SQL and R) but packages/tools do not
>> always adhere to this...so it's never clear how the output will be
>> ordered.  Sorry for the confusion.
>>
>>
>> On Tue, Nov 26, 2019 at 12:22 PM Leonidas Liakos <
>> [hidden email]> wrote:
>>
>>> Why do they seem logical since they do not match?
>>>
>>> Check for example index 1 (Sunday). The results are different for the
>>> three processes
>>>
>>> > stackapply_mean
>>> class      : RasterBrick
>>> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>> resolution : 500, 500  (x, y)
>>> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>> crs        : NA
>>> source     :
>>> /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
>>> names      :  index_5,  index_6,  index_7,  index_1,  index_2,
>>> index_3,  index_4
>>> min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
>>> 442.7436
>>> max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
>>> 556.2471
>>>
>>>
>>> > ver_mean
>>> class      : RasterStack
>>> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>> resolution : 500, 500  (x, y)
>>> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>> crs        : NA
>>> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,
>>> layer.6,  layer.7
>>> min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
>>> 429.6900
>>> max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
>>> 561.7972
>>>
>>>
>>> > z
>>> class      : RasterBrick
>>> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>> resolution : 500, 500  (x, y)
>>> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>> crs        : NA
>>> source     :
>>> /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
>>> names      :       X1,       X2,       X3,       X4,       X5,
>>> X6,       X7
>>> min values : 440.0467, 444.9182, 437.1589, 444.6946, 440.2028, 429.6900,
>>> 442.7436
>>> max values : 563.8341, 561.7687, 560.4509, 565.8671, 560.1375, 561.7972,
>>> 556.2471
>>>            : 1, 2, 3, 4, 5, 6, 7
>>>
>>>
>>> On 11/26/19 7:03 PM, Vijay Lulla wrote:
>>>
>>> If you read the code/help for `stackApply` and `zApply` you'll see that
>>> the results that you obtain make sense (at least they seem
>>> sensible/reasonable to me).  IMO, if you want to control the ordering of
>>> your layers then just use sapply, like how you've used for ver_mean.  IMO,
>>> this is the only reliable (safe?), and quite a readable, way to accomplish
>>> what you're trying to do.
>>> Just my 2 cents.
>>> -- Vijay.
>>>
>>> On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via R-sig-Geo <
>>> [hidden email]> wrote:
>>>
>>>> I added raster::zApply in my tests to validate the results. However, the
>>>> indices of the names of the results are different now. Recall that the
>>>> goal is to calculate from a raster stack time series the mean per day of
>>>> the week. And that problem I have is that stackApply, zApply and
>>>> calc/sapply return different indices in the result names. New code is
>>>> available here:
>>>> https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>>>> I'm really curious about missing something.
>>>>
>>>>
>>>> On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>>>> > 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 4:15 PM Leonidas Liakos via R-sig-Geo
>>>> > <[hidden email] <mailto:[hidden email]>> wrote:
>>>> >
>>>> >     I run the example with clusterR:
>>>> >
>>>> >     no_cores <- parallel::detectCores() -1
>>>> >     raster::beginCluster(no_cores)
>>>> >     ?????? res <- raster::clusterR(inp, raster::stackApply, args =
>>>> >     list(indices=c(2,2,3,3,1,1),fun = mean))
>>>> >     raster::endCluster()
>>>> >
>>>> >     And the result is:
>>>> >
>>>> >     > res
>>>> >     class?????????? : RasterBrick
>>>> >     dimensions : 180, 360, 64800, 3?? (nrow, ncol, ncell, nlayers)
>>>> >     resolution : 1, 1?? (x, y)
>>>> >     extent???????? : -180, 180, -90, 90?? (xmin, xmax, ymin, ymax)
>>>> >     crs?????????????? : +proj=longlat +datum=WGS84 +ellps=WGS84
>>>> >     +towgs84=0,0,0
>>>> >     source???????? : memory
>>>> >     names?????????? : layer.1, layer.2, layer.3
>>>> >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>>>> >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>>> >
>>>> >
>>>> >     layer.1, layer.2, layer.3 (?)
>>>> >
>>>> >     So what corrensponds to what?
>>>> >
>>>> >
>>>> >     If I run:
>>>> >
>>>> >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>>> >
>>>> >     The result is:
>>>> >
>>>> >     > res2
>>>> >     class      : RasterBrick
>>>> >     dimensions : 180, 360, 64800, 3  (nrow, ncol, ncell, nlayers)
>>>> >     resolution : 1, 1  (x, y)
>>>> >     extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>>>> >     crs        : +proj=longlat +datum=WGS84 +ellps=WGS84
>>>> +towgs84=0,0,0
>>>> >     source     : memory
>>>> >     names      : index_2, index_3, index_1
>>>> >     min values :     1.5,     3.5,     5.5
>>>> >     max values :     1.5,     3.5,     5.5
>>>> >
>>>> >     There is no consistency with the names of the output and obscure
>>>> >     correspondence with the indices in the case of clusterR
>>>> >
>>>> >
>>>> >             [[alternative HTML version deleted]]
>>>> >
>>>> >     _______________________________________________
>>>> >     R-sig-Geo mailing list
>>>> >     [hidden email] <mailto:[hidden email]>
>>>> >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> >
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>
>>>
>>
>

--
Vijay Lulla

Assistant Professor,
Dept. of Geography, IUPUI
425 University Blvd, CA-207C.
Indianapolis, IN-46202
[hidden email]
ORCID: https://orcid.org/0000-0002-0823-2522
Online: http://vijaylulla.com | https://github.com/vlulla

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: raster: stackApply problems..

Roger Bivand
Administrator
In reply to this post by R-sig-geo mailing list
On Wed, 27 Nov 2019, Leonidas Liakos via R-sig-Geo wrote:

> Thank you for your help!
>
> I tried to fix stackApply according to your instructions.
>
> Now the indices of names are the same and consistent with indices
> enumeration (gist for validation and tests:
> https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170).
>
> I've attached a patch file here:
>
> https://gist.github.com/kokkytos/ca2c319134677b19900579665267a7a7
Thanks very much for contributing!

Please consider raising an issue on https://github.com/rspatial/raster 
pointing to this thread and your patch. I had expected response from
raster developers here, but they may well be on field work, so raising an
issue on the development site should get their attention when there is
enough time for them to look. You might even fork raster, apply your patch
and file a PR in addition to the issue. In that case, a short test would
be helpful, and maybe edits to the documentation.

Roger


>
>> stackapply_mean
> class      : RasterBrick
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> source     : /tmp/Rtmp9W8UNc/raster/r_tmp_2019-11-27_191205_2929_20324.grd
> names      :  index_5,  index_6,  index_7,  index_1,  index_2, 
> index_3,  index_4
> min values : 444.6946, 440.2028, 429.6900, 442.7436, 440.0467, 444.9182,
> 437.1589
> max values : 565.8671, 560.1375, 561.7972, 556.2471, 563.8341, 561.7687,
> 560.4509
>
>> ver_mean
> class      : RasterStack
> dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
> resolution : 500, 500  (x, y)
> extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
> crs        : NA
> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5, 
> layer.6,  layer.7
> min values : 442.7436, 440.0467, 444.9182, 437.1589, 444.6946, 440.2028,
> 429.6900
> max values : 556.2471, 563.8341, 561.7687, 560.4509, 565.8671, 560.1375,
> 561.7972
>
>
>
>
>
> On 11/26/19 10:58 PM, Vijay Lulla wrote:
>> Hmm...it appears that stackApply is using different conditions which
>> might be causing this problem. Below is the snippet of the code which
>> I think might be the problem.
>>
>> ## For canProcessInMemory
>> if (rowcalc) {
>>   v <- lapply(uin, function(i) fun(x[, ind == i, drop = FALSE], na.rm
>> = na.rm))
>> }
>> else {
>>   v <- lapply(uin, function(i, ...) apply(x[, ind == i, drop = FALSE],
>> 1, fun, na.rm = na.rm))
>> }
>>
>>
>> ## If canProcessInMemory is not TRUE
>> if (rowcalc) {
>>   v <- lapply(uin, function(i) fun(a[, ind == uin[i], drop = FALSE],
>> na.rm = na.rm))
>> }
>> else {
>>   v <- lapply(uin, function(i, ...) apply(a[, ind == uin[i], drop =
>> FALSE], 1, fun, na.rm = na.rm))
>> }
>>
>> I think they should both be same but it appears that they aren't and
>> that's what you've discovered.  Maybe you can try fix(stackApply) to
>> see if it really is the problem and can tell us what you find. 
>> Anyways, good catch...and...sorry for wasting your time.
>> Cordially,
>> Vijay.
>>
>> On Tue, Nov 26, 2019 at 2:53 PM Leonidas Liakos
>> <[hidden email] <mailto:[hidden email]>> wrote:
>>
>>     Thank you!
>>     The problem is not with the resulting values but with the index
>>     mapping. Values are correct in all three cases.
>>
>>     As I wrote in a previous post in the thread
>>     (https://stat.ethz.ch/pipermail/r-sig-geo/2019-November/027821.html)
>>     , stackApply behaves inconsistently depending on whether the
>>     exported stack will remain in memory or it will be stored, due to
>>     its large size, on the hard disk.
>>
>>     In the first case the indices are identical to my function
>>     (ver_mean) and the lubridate::wday indexing system (and are
>>     correct) while in the second they are shuffled.
>>
>>     So, while Sunday has index 1 and while in the first case (when the
>>     result is in memory) stackApply returns the correct index, in the
>>     second case (when the result is stored on the hard disk) it
>>     returns index_4! So how can one be sure if index_1 corresponds to
>>     Sunday or another day using stackApply since it sometimes
>>     enumerates it with index_1 and sometimes index_4?
>>
>>
>>     Try to run this example (when the resulting stack remains in
>>     memory) to see that the indexes are identical (stackApply =
>>     ver_median = lubridate::wday)
>>     https://gist.github.com/kokkytos/5d554b5a725bb48d2189e2d1fa0e2206
>>
>>     Thank you again
>>
>>     On 11/26/19 9:00 PM, Vijay Lulla wrote:
>>>     I'm sorry for the miscommunication.  What I meant to say is that
>>>     the output from stackApply and zApply are the same (because
>>>     zApply uses stackApply internally) except the names.  The
>>>     behavior of stackApply makes sense because AFAIUI R doesn't
>>>     automatically reorder vectors/indices that it receives.  Your
>>>     observation about inconsistent result with ver_mean is very valid
>>>     though!  And, that's what I meant with my comment that using
>>>     sapply with the explicit ordering that you want is the best way
>>>     to control what raster package will output.  In R the input order
>>>     should be maintained (this is the prime difference between SQL
>>>     and R) but packages/tools do not always adhere to this...so it's
>>>     never clear how the output will be ordered.  Sorry for the confusion.
>>>
>>>
>>>     On Tue, Nov 26, 2019 at 12:22 PM Leonidas Liakos
>>>     <[hidden email] <mailto:[hidden email]>> wrote:
>>>
>>>         Why do they seem logical since they do not match?
>>>
>>>         Check for example index 1 (Sunday). The results are different
>>>         for the three processes
>>>
>>>        > stackapply_mean
>>>         class      : RasterBrick
>>>         dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>>         resolution : 500, 500  (x, y)
>>>         extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>>         crs        : NA
>>>         source     :
>>>         /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191359_7710_20324.grd
>>>         names      :  index_5,  index_6,  index_7,  index_1, 
>>>         index_2,  index_3,  index_4
>>>         min values : 440.0467, 444.9182, 437.1589, 444.6946,
>>>         440.2028, 429.6900, 442.7436
>>>         max values : 563.8341, 561.7687, 560.4509, 565.8671,
>>>         560.1375, 561.7972, 556.2471
>>>
>>>
>>>        > ver_mean
>>>         class      : RasterStack
>>>         dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>>         resolution : 500, 500  (x, y)
>>>         extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>>         crs        : NA
>>>         names      :  layer.1,  layer.2,  layer.3,  layer.4, 
>>>         layer.5,  layer.6,  layer.7
>>>         min values : 442.7436, 440.0467, 444.9182, 437.1589,
>>>         444.6946, 440.2028, 429.6900
>>>         max values : 556.2471, 563.8341, 561.7687, 560.4509,
>>>         565.8671, 560.1375, 561.7972
>>>
>>>
>>>        > z
>>>         class      : RasterBrick
>>>         dimensions : 300, 300, 90000, 7  (nrow, ncol, ncell, nlayers)
>>>         resolution : 500, 500  (x, y)
>>>         extent     : 0, 150000, 0, 150000  (xmin, xmax, ymin, ymax)
>>>         crs        : NA
>>>         source     :
>>>         /tmp/RtmpkRMXLb/raster/r_tmp_2019-11-26_191439_7710_04780.grd
>>>         names      :       X1,       X2,       X3,       X4,      
>>>         X5,       X6,       X7
>>>         min values : 440.0467, 444.9182, 437.1589, 444.6946,
>>>         440.2028, 429.6900, 442.7436
>>>         max values : 563.8341, 561.7687, 560.4509, 565.8671,
>>>         560.1375, 561.7972, 556.2471
>>>                    : 1, 2, 3, 4, 5, 6, 7
>>>
>>>
>>>         On 11/26/19 7:03 PM, Vijay Lulla wrote:
>>>>         If you read the code/help for `stackApply` and `zApply`
>>>>         you'll see that the results that you obtain make sense (at
>>>>         least they seem sensible/reasonable to me).  IMO, if you
>>>>         want to control the ordering of your layers then just use
>>>>         sapply, like how you've used for ver_mean.  IMO, this is the
>>>>         only reliable (safe?), and quite a readable, way to
>>>>         accomplish what you're trying to do.
>>>>         Just my 2 cents.
>>>>         -- Vijay.
>>>>
>>>>         On Tue, Nov 26, 2019 at 11:19 AM Leonidas Liakos via
>>>>         R-sig-Geo <[hidden email]
>>>>         <mailto:[hidden email]>> wrote:
>>>>
>>>>             I added raster::zApply in my tests to validate the
>>>>             results. However, the
>>>>             indices of the names of the results are different now.
>>>>             Recall that the
>>>>             goal is to calculate from a raster stack time series the
>>>>             mean per day of
>>>>             the week. And that problem I have is that stackApply,
>>>>             zApply and
>>>>             calc/sapply return different indices in the result
>>>>             names. New code is
>>>>             available here:
>>>>             https://gist.github.com/kokkytos/93f315a5ecf59c0b183f9788754bc170
>>>>             I'm really curious about missing something.
>>>>
>>>>
>>>>             On 11/20/19 3:30 AM, Frederico Faleiro wrote:
>>>>            > 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 4:15 PM Leonidas Liakos via
>>>>             R-sig-Geo
>>>>            > <[hidden email]
>>>>             <mailto:[hidden email]>
>>>>             <mailto:[hidden email]
>>>>             <mailto:[hidden email]>>> wrote:
>>>>            >
>>>>            >     I run the example with clusterR:
>>>>            >
>>>>            >     no_cores <- parallel::detectCores() -1
>>>>            >     raster::beginCluster(no_cores)
>>>>            >     ?????? res <- raster::clusterR(inp,
>>>>             raster::stackApply, args =
>>>>            >     list(indices=c(2,2,3,3,1,1),fun = mean))
>>>>            >     raster::endCluster()
>>>>            >
>>>>            >     And the result is:
>>>>            >
>>>>            >     > res
>>>>            >     class?????????? : RasterBrick
>>>>            >     dimensions : 180, 360, 64800, 3?? (nrow, ncol,
>>>>             ncell, nlayers)
>>>>            >     resolution : 1, 1?? (x, y)
>>>>            >     extent???????? : -180, 180, -90, 90?? (xmin, xmax,
>>>>             ymin, ymax)
>>>>            >     crs?????????????? : +proj=longlat +datum=WGS84
>>>>             +ellps=WGS84
>>>>            >     +towgs84=0,0,0
>>>>            >     source???????? : memory
>>>>            >     names?????????? : layer.1, layer.2, layer.3
>>>>            >     min values :???????? 1.5,???????? 3.5,???????? 5.5
>>>>            >     max values :???????? 1.5,???????? 3.5,???????? 5.5??
>>>>            >
>>>>            >
>>>>            >     layer.1, layer.2, layer.3 (?)
>>>>            >
>>>>            >     So what corrensponds to what?
>>>>            >
>>>>            >
>>>>            >     If I run:
>>>>            >
>>>>            >     res2 <- stackApply(inp,c(2,2,3,3,1,1),mean)
>>>>            >
>>>>            >     The result is:
>>>>            >
>>>>            >     > res2
>>>>            >     class      : RasterBrick
>>>>            >     dimensions : 180, 360, 64800, 3  (nrow, ncol,
>>>>             ncell, nlayers)
>>>>            >     resolution : 1, 1  (x, y)
>>>>            >     extent     : -180, 180, -90, 90  (xmin, xmax,
>>>>             ymin, ymax)
>>>>            >     crs        : +proj=longlat +datum=WGS84
>>>>             +ellps=WGS84 +towgs84=0,0,0
>>>>            >     source     : memory
>>>>            >     names      : index_2, index_3, index_1
>>>>            >     min values :     1.5,     3.5,     5.5
>>>>            >     max values :     1.5,     3.5,     5.5
>>>>            >
>>>>            >     There is no consistency with the names of the
>>>>             output and obscure
>>>>            >     correspondence with the indices in the case of
>>>>             clusterR
>>>>            >
>>>>            >
>>>>            >             [[alternative HTML version deleted]]
>>>>            >
>>>>            >     _______________________________________________
>>>>            >     R-sig-Geo mailing list
>>>>            >     [hidden email]
>>>>             <mailto:[hidden email]>
>>>>             <mailto:[hidden email]
>>>>             <mailto:[hidden email]>>
>>>>            >     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>            >
>>>>
>>>>                     [[alternative HTML version deleted]]
>>>>
>>>>             _______________________________________________
>>>>             R-sig-Geo mailing list
>>>>             [hidden email] <mailto:[hidden email]>
>>>>             https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>>
>>>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway