Re: [R-sig-Geo] raster: stackApply problems..
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, jon.sko...@ec.europa.eu wrote:
>
> Leonidas,
>
> I see that you are not happy with the output, but it is not so clear
> what you actually expect to see.
>
>
> If you use stackApply directly, the indices are used in the names.
> Layer 1 and 8 belong to the group with index 4. It is the first group
> in the list of indexes, so the first layer of the output is then
> referred to as index_4. Then comes index_5 with layers 2, 10 and 15 of
> your input. The order of these names will follow the order of the
> first appearance of your indices. The indices gets lost with the use
> of clusterR, so it gives you the same output, but with names layer.1 -
> layer.7.
>
>
> You could change the names of the result from clusterR with:
>
> names(ResClusterR) = paste0("index_", unique(indices))
>
>
> If you want your result (from stackApply or clusterR) to follow the
> order of your indices, you should be able to get this with:
>
>
> sResClusterR = ResClusterR[[order(names(ResClusterR))]]
>
>
> Does this help you further?
>
> Jon
>
>
>
>
> --
> Jon Olav Skøien
> European Commission
> Joint Research Centre – JRC.E.1
> Disaster Risk Management Unit
> Building 26b 1/144 | Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 267
> jon.sko...@ec.europa.eu Tel: +39 0332 789205 Disclaimer: Views
> expressed in this email are those of the individual and do not
> necessarily represent official views of the European Commission.
>
>
>
>
> *From:* R-sig-Geo on behalf of
> Leonidas Liakos via R-sig-Geo
> *Sent:* 21 November 2019 08:52
> *To:* Ben Tupper; r-sig-geo@r-project.org
> *Subject:* Re: [R-sig-Geo] raster: stackApply problems..
>
> Unfortunately the names are not always in ascending order. This is the
> result of my data.
>
> names : index_4, index_5, index_6, index_7, index_1, index_2, index_3
> min values : 3, 3, 3, 3, 3, 3, 3
> max values : 307.0, 297.5, 311.0, 313.0, 468.0, 290.0, 302.0
>
> And worst of all, it is not a proper match with indices.
>
> If I run it with clusterR then the result is different:
>
> names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7
> min values : 3, 3, 3, 3, 3, 3, 3
> max values : 307.0, 297.5, 311.0, 313.0, 468.0, 290.0,
> 302.0
>
>
> The solution is to reorder the layers of the stack so that the
> stackApply indices are in ascending order e.g. 1,1,1,2,2,2,3,3,3 ...
>
> My indices of my data was like that:
>
> 4 5 6 7 1 2 3 4 5 6 7 1 2 3 5 6 7
>
> I've reported this behavior here
> https://urldefense.com/v3/__https://github.com/rspatial/raster/issues/82__;!NW73rmyV52c!SiZfwLn8F-IC_xeeUNNjzf8STJX1LMbYaoJKqfWo5ImGWi_dEhB7ilEG9835G-KjNZJC$
>
>
>
> On 11/20/19 3:05 PM, Ben Tupper wrote:
> > Hi,
> >
> > That is certainly is unexpected to have two different naming styles.
> > It's not really solution to take to the bank, but you could simply
> > compose your own names assuming that the layer orders are always
> > returned in ascending index order.
> > Would that work for you
> >
> > ### start
> > library(raster)
> >
> > # Compute layer names for stackApply output
> > #
> > # @param index numeric, 1-based layer indices used for stackApply
> function
> > # @param prefix character, prefix for names
> > # @return character layers names in index order
> > layer_names <- function(index = c(2,2,3,3,1,1), prefix = c("layer.",
> > "index_")[1]){
> > paste0(prefix, sort(unique(index)))
> > }
> >
> > indices <- c(2,2,3,3,1,1)
> >
> > r <- raster()
> > values(r) <- 1
> > # simple sequential stack from 1 to 6 in all cells
> > s <- stack(r, r*2, r*3, r*4, r*5, r*6)
> > s
> >
> > beginCluster(2)
> > res <- clusterR(s, stackApply, args = list(indices=indices, fun = mean))
> > raster::endCluster()
> > names(res) <- layer_names(indices, prefix = "foobar.")
> > res
> >
> > res2 <- stackApply(s, indices, mean)
> > names(res2) <- layer_names(indices, prefix = "foobar.")
> > res2
> > ### end
> >
> >
> > On Wed, Nov 20, 2019 at 1:36 AM Leonidas Liakos via R-sig-Geo
> > wrote:
> >> This is not a reasonable solution. It is not efficient to run
> stackapply
> >> twice to get
Re: [R-sig-Geo] sp/rgdal workflows with PROJ >= 6 and GDAL >= 3
A description of the status now with regard to a prototype resolution is online at: https://rsbivand.github.io/ECS530_h19/ECS530_III.html I'm planning to stream a talk about this at 09:15-11:00 CET on Tuesday 3 December. I need a volunteer to test the streaming link in advance during next week. I'm unsure which technology to use for remote participants to provide feedback. Contributions/comments welcome! Roger On Fri, 15 Nov 2019, Roger Bivand wrote: The development version of rgdal on R-Forge is now at rev 894, and is now ready for trying out with PROJ6/GDAL3 workflows, and workflows that may migrate within 6 months to modern CRS representations. The motivating RFC is also updated to cover coordinate operations, the use of prepared (pre-searched) coordinate operations, and should be read carefully by anyone using rgdal::spTransform(). Note further that rgdal::project() will not be adapted for PROJ6, and is effectively deprecated. I'll be running reverse dependency checks, and may be bugging package maintainers. I would really prefer that mainainers of packages using spTransform() checked themselves and joined this thread or the associated twitter thread: https://twitter.com/RogerBivand/status/1194586193108914177 Be ready for modern PROJ and GDAL, they are already being deployed across open source geospatial software, like GRASS, QGIS, pyproj, spatialite etc. Waiting, hopefully not in vain, for contributions. Roger On Wed, 13 Nov 2019, Roger Bivand wrote: And this link explains the CDN proposal for grid distribution: https://www.spatialys.com/en/crowdfunding/ Roger On Wed, 13 Nov 2019, Roger Bivand wrote: Because PROJ >= 6 and GDAL >= 3 change the way that PROJ strings (representations of coordinate reference systems) are handled, steps are being taken to find ways to adapt sp/rgdal workflows. A current proposal is to store the WKT2_2018 string as a comment to CRS objects as defined in the sp package. A draft development-in-progress version of rgdal is available at https://r-forge.r-project.org/R/?group_id=884, and for sp at https://github.com/rsbivand/sp (this version of sp requires rgdal >= 1.5-1). This adds the WKT comments to CRS objects on reading vector and raster data sources, and uses WKT comments if found when writing vector and raster objects (or at least does as far as I've checked, possibly fragile). An RFC with tersely worked cases for using CRS object comments to carry WKT strings but maintaining full backward compatibility is online at http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html. If you have other ideas or concerns about trying to use this mechanism for sp CRS objects, please contribute at your earliest convenience. http://rgdal.r-forge.r-project.org/reference/list_coordOps.html shows the beginning of the next step, to query transformation operations to find viable coordinate operation pipelines. I'm assuming that the previous behaviour (transform without considering accuracy with whatever is to hand) is not viable going forward, and that we will need two steps: list coordinate operations between source and target CRS (using the WKT comments as better specifications than the PROJ strings), possibly intervene manually to install missing grids, then undertake the coordinate operation. The fallback may be simply to choose the least inaccurate available coordinate operation, but this should be a fallback. This means that all uses of spTransform() will require intervention. Is this OK (it is tiresome but modernises workflows once), or is it not OK (no user intervention is crucial)? These behaviours may be set in an option, so that package maintainers and users may delay modernisation, but all are undoubtedly served by rapid adaptation (GRASS 7.8.1 released yesterday, libspatialite, pyproj, QGIS development versions all state that they list candidate coordinate operations). We cannot ship all the grids, they are very bulky, and probably nobody needs sub-metre accuracy world-wide. Work in PROJ is starting to create a content delivery network for trusted download and mechanisms for registering downloaded grids on user platforms. We would for example not want Windows users of rgdal and sf to have to download the same grid twice. Comments welcome here and at https://github.com/r-spatial/discuss/issues/28 or https://github.com/r-spatial/sf/issues/1187 Roger -- Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no https://orcid.org/-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0J=en ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
