[R-sig-Geo] raster: stackApply problems..
Leonidas Liakos
|eon|d@@_||@ko@ @end|ng |rom y@hoo@gr
Fri Nov 22 20:53:12 CET 2019
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, Jon.SKOIEN using ec.europa.eu wrote:
>
> Leonidas,
>
> I see that you are not happy with the output, but it is not so clear
> what you actually expect to see.
>
>
> If you use stackApply directly, the indices are used in the names.
> Layer 1 and 8 belong to the group with index 4. It is the first group
> in the list of indexes, so the first layer of the output is then
> referred to as index_4. Then comes index_5 with layers 2, 10 and 15 of
> your input. The order of these names will follow the order of the
> first appearance of your indices. The indices gets lost with the use
> of clusterR, so it gives you the same output, but with names layer.1 -
> layer.7.
>
>
> You could change the names of the result from clusterR with:
>
> names(ResClusterR) = paste0("index_", unique(indices))
>
>
> If you want your result (from stackApply or clusterR) to follow the
> order of your indices, you should be able to get this with:
>
>
> sResClusterR = ResClusterR[[order(names(ResClusterR))]]
>
>
> Does this help you further?
>
> Jon
>
>
>
>
> --
> Jon Olav Skøien
> European Commission
> Joint Research Centre – JRC.E.1
> Disaster Risk Management Unit
> Building 26b 1/144 | Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 267
> jon.skoien using ec.europa.eu Tel: +39 0332 789205 Disclaimer: Views
> expressed in this email are those of the individual and do not
> necessarily represent official views of the European Commission.
>
>
>
> ------------------------------------------------------------------------
> *From:* R-sig-Geo <r-sig-geo-bounces using r-project.org> on behalf of
> Leonidas Liakos via R-sig-Geo <r-sig-geo using r-project.org>
> *Sent:* 21 November 2019 08:52
> *To:* Ben Tupper; r-sig-geo using r-project.org
> *Subject:* Re: [R-sig-Geo] raster: stackApply problems..
>
> Unfortunately the names are not always in ascending order. This is the
> result of my data.
>
> names : index_4, index_5, index_6, index_7, index_1, index_2, index_3
> min values : 3, 3, 3, 3, 3, 3, 3
> max values : 307.0, 297.5, 311.0, 313.0, 468.0, 290.0, 302.0
>
> And worst of all, it is not a proper match with indices.
>
> If I run it with clusterR then the result is different:
>
> names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7
> min values : 3, 3, 3, 3, 3, 3, 3
> max values : 307.0, 297.5, 311.0, 313.0, 468.0, 290.0,
> 302.0
>
>
> The solution is to reorder the layers of the stack so that the
> stackApply indices are in ascending order e.g. 1,1,1,2,2,2,3,3,3 ...
>
> My indices of my data was like that:
>
> 4 5 6 7 1 2 3 4 5 6 7 1 2 3 5 6 7
>
> I've reported this behavior here
> https://urldefense.com/v3/__https://github.com/rspatial/raster/issues/82__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G-KjNZJC$
>
>
>
> On 11/20/19 3:05 PM, Ben Tupper wrote:
> > Hi,
> >
> > That is certainly is unexpected to have two different naming styles.
> > It's not really solution to take to the bank, but you could simply
> > compose your own names assuming that the layer orders are always
> > returned in ascending index order.
> > Would that work for you
> >
> > ### start
> > library(raster)
> >
> > # Compute layer names for stackApply output
> > #
> > # @param index numeric, 1-based layer indices used for stackApply
> function
> > # @param prefix character, prefix for names
> > # @return character layers names in index order
> > layer_names <- function(index = c(2,2,3,3,1,1), prefix = c("layer.",
> > "index_")[1]){
> > paste0(prefix, sort(unique(index)))
> > }
> >
> > indices <- c(2,2,3,3,1,1)
> >
> > r <- raster()
> > values(r) <- 1
> > # simple sequential stack from 1 to 6 in all cells
> > s <- stack(r, r*2, r*3, r*4, r*5, r*6)
> > s
> >
> > beginCluster(2)
> > res <- clusterR(s, stackApply, args = list(indices=indices, fun = mean))
> > raster::endCluster()
> > names(res) <- layer_names(indices, prefix = "foobar.")
> > res
> >
> > res2 <- stackApply(s, indices, mean)
> > names(res2) <- layer_names(indices, prefix = "foobar.")
> > res2
> > ### end
> >
> >
> > On Wed, Nov 20, 2019 at 1:36 AM Leonidas Liakos via R-sig-Geo
> > <r-sig-geo using r-project.org> 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 <
> >>> r-sig-geo using r-project.org> 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
> >>>> R-sig-Geo using r-project.org
> >>>>
> 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
> >> R-sig-Geo using r-project.org
> >>
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
> >
> >
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-geo__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G3J81kzN$
>
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list