is it possible to support tile projections other than 4326?

[jmlondon]

My reading of the documentation suggests that, while non-4326 rasters are supported, they are transformed to 4326 as part of the tile function. So, there is not currently a way to generate tiles in other projections.

Since the leaflet package now supports projections other than Web Mercator (via the Proj4Leaflet add-in), it would be beneficial if tiler could create tiles in other projections.

This is of particular interest for maps with data that cross the anti-meridian (180 longitude) or in polar regions which leaftlet (Web Mercator) currently struggles with.

[leonawicz]

Totally, I want to do that. I don't like having just the two options when really any projection should be possible and Leaflet supports them all. I chose the current capabilities as the initial low-hanging fruit. I'm guessing it's (hopefully) as simple as telling gdal2tiles via command line argument inside tile what the output projection should be, but I haven't peered that closely or had time to really think about it yet. If you happen to already know this please let me know. I make a lot of maps of Alaska and really dislike the Web Mercator view even aside from the anti-meridian issue.

I remember playing with other projections in leaflet using tutorial examples I found online (this was before I made tiler) and they worked great, but for whatever reason my attempts to use maps in NAD83 Alaska Albers EAC projection always came out looking really awful and I couldn't figure out what I was doing wrong, if anything. So at the time I didn't feel very confident about it.

But yeah, ideally it will be as simple as not forcing a reprojection unless the desired output CRS is different from that of the input file. Also, I could be wrong on this, but I have a nagging feeling that gdal2tiles might just want EPSG codes, not proj4 strings like are attached to a raster object.

[jmlondon]

I managed to get the ESRI Arctic Ocean Basemap working in leaflet a few weeks ago (see rstudio/leaflet#550). If you're mapping north of the Aleutians, it's a pretty good option. The code is a bit clunky so I hope to find some time and

I have no experience with gdal2tiles except I looked into it a while back when I had an idea to re-project the natural earth map tiles (naturalearthdata.com) into something more Alaska/N. Pacific friendly. I dropped it b/c I realized I didn't have a way to serve up the tiles, but I like your idea of pushing to github. Maybe I'll revive that idea this week.

[leonawicz]

I checked and gdal2tiles has an SRS argument, which I remembered seemed to pertain to input, not output. I did some testing and, so far as I can tell, it continues to appear to be an unnecessary argument, as gdal2tiles is already picking up the CRS of the input raster.

Using a new local repo branch I've tested an example using a projected Alaska map. Providing SRS internally to gdal2tiles or continuing to leave it out had no apparent effect. However, by simply no longer forcing the reprojection to EPSG 4326, I did get a different map out. So, it may be that easy, or close to that easy.

The problem is that the tiles appear to be the opposite of what I'd expect. Instead of warping Alaska, it still appears exactly as before, looking like the awkward WGS84 Alaska map, but instead the background of the map domain is curved. Maybe I am telling it to do something backwards, but I don't know at this point.

[jmlondon]

I think the default profile for gdal2tiles is mercator. So, maybe, specify --profile=raster ??

I'm playing around with things for a bit today as well. I'll report back if I have any luck.

[jmlondon]

I got this via the following:

gdal_translate -a_srs EPSG:4326 -gcp 0 0 -180 90 -gcp 16200 0 180 90 -gcp 16200 8100 180 -90 NE1_HR_LC_SR_W.tif NE1_HR_LC_SR_W_proj.tif
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3571 -co TILED=YES NE1_HR_LC_SR_W.tif naturalearth_3571.tif
gdal2tiles.py --profile=raster -z 1-9 naturalearth_3571.tif naturalearth_3571

Obviously, the tearing at 180 isn't cool. But, I'm hoping there's a better solution available. The key is, I think the --profile=raster option to gda2tiles.py resulted in the creation of tiles in the 3571 projection.

[leonawicz]

That was totally my issue! Heh, I even overlooked that I had -p raster in the internal call to the alternative gdal2tiles script for the non-geographic case (these eyes get tired).

My quick example is not as robust as yours. I can't tell where the 180 is, but it's looking pretty good for what it is. I think it's actually all there. I can push this projections test branch to GitHub soon if you want to try it out.

Definitely need the default Mercator for the usual Leaflet use cases (like reprojecting this map to 4326), so I'm looking at using a new argument, profile to pass along to gdal2tiles. At the moment, I'm thinking the least annoying combination of defaults is:

(1) no reprojection of file before tiling unless the user specifies an output CRS with crs_out (crs_out = NULL the default). This keeps a perspective in tiler of not assuming what the user wants. Unless otherwise specified, tiles are made based on the input file (or crs override) without altering things such as assuming the tiles are for Google Maps or something.

(2) mercator is still the default tile cutting profile, just like in gdal2tiles, and there's no reason to try and guess what would be the most likely CRS of input files people would provide. If file is not in a Mercator CRS, then the user can add profile = "raster" for example.

Even if EPSG 4326 is the dominant use case output, regardless of input, I'm leaning toward this perspective of defaulting to deriving output CRS from input CRS and still requiring users to specify crs_out for EPSG 4326. Since tiler is about making tiles, not about making a specific type of popular map layout.

Let me know if you have thoughts on this. It's hard to guess how everyone else will see it. But I think this is a good perspective to retain for the package. The only noticeable "big" change from the current CRAN version would be no more forced reprojection to 4326. Reprojection becoming an option via crs_out, with a little help from profile as needed.

[leonawicz]

Okay, try out this one: devtools::install_github("leonawicz/tiler", ref = "projections")

The projections branch makes the changes I mentioned. I've kept everything working technically, but that doesn't mean for sure that different projections are all working correctly. This is looking promising though.

I also hate lacking documentation, so in this branch, the help docs, vignette and NEWS files are updated to match the code changes I made to tile. So feel free to refer to them as needed regarding the new tile arguments. But overall I really didn't change much. Let me know what you find. Thanks! :)

[jmlondon]

I'll give this a try later today or tomorrow. For completeness, I think I can fix my gdalwarp issue with the natural earth dataset by specifying -wo SOURCE_EXTRA=1000. This essentially creates some overlap and prevents the ugly tearing at 180.

Once I have that finalized, I'll give tiler a try to actually make the tiles and then host on my github for inclusion in calls to the leaflet package. For this, I'll need to sort out the specifications for leafletCRS(), but I'm hoping that won't be too hard.

[jmlondon]

I was able to install the "projections" branch and gave it a try with the following code

library(tiler)
map <- "~/Downloads/NE2_HR_LC_SR_W_DR/naturalearth_3571.tif"
tile_dir <- "tiles"
tile(map, tile_dir, "5-7", profile = "raster")

That resulted in the following warnings/error

partial argument match of 'along' to 'along.with'
partial argument match of 'along' to 'along.with'
partial argument match of 'along' to 'along.with'
partial argument match of 'along' to 'along.with'
partial argument match of 'along' to 'along.with'
partial argument match of 'along' to 'along.with'
partial argument match of 'along' to 'along.with'
Preparing for tiling...
Error: vector memory exhausted (limit reached?)

I am able to successfully create tiles directly with gdal2tiles.py. So, I'm not sure why vector memory was exhausted. I'm also not sure what the partial argument match ... error is all about.

[leonawicz]

I thought I just downloaded the same file, but when I unzipped it it turned out to be named NE2_HR_LC_SR_W_DR.tif. It's also a nearly 700 MB RGB raster. Are you using something different? I'd like to reproduce what you attempted. I haven't tried tiling anything near this large.

[jmlondon]

I changed things up a bit ...

I have been working with the NE2_HR_LC_SR_W_DR.zip which is the large Natural Earth II dataset.

It certainly takes time (~20-30 minutes) to work through the process. But, most of that is in the gdalwarp stage and because I was going for the lanczos option.

I'll update my repository on github (probably later this evening) and, hopefully, it will provide some insight so you can try and reproduce things on your end.

[leonawicz]

I pushed a tiny update (same projections branch) that seems to help with memory management on the R portion of the process. At least, it did on my Windows machine, which was otherwise maxing out my RAM with such a big file (and going to disk, so not failing). I did the tiling on zoom 0-7 and it worked great either way, just that this new way things looked a lot more stable.

This was only using the default CRS of the input file. I did not reproject anything with crs_out.

My recommendation in general, but which I guess would apply even more to really large files, is that file should ideally be prepped to represent exactly what is to be tiled. I.e., reproject it first, or even make a RGB raster out of it first if it's not already RGB-layered. Basically, anything that the raster package is going to have to do to file inside the tile call is going to bring that data into memory in R and then write the result to a temp file and the temp file will be read by gdal2tiles for tiling (and then the temp file is deleted). So if the user really wants to avoid overhead there with a really big file, this is the better approach in such a situation. If file doesn't need any pre-processing, R will never have to load it, and file itself can be passed to gdal2tiles instead of an intermediary file created by R.

Let me know if you suggest trying out a particular Proj4 string for crs_out.

[jmlondon]

My first crack at documenting my process is available here.

I do exactly what you said and take care of all the re-projecting via gdalwarp before passing that to tiler. Note, the current repo doesn't include use of your updated fix. So, I'll give that a try.

Also, note that my attempts to get the tile set working with leaflet is also currently failing.

[leonawicz]

I was able to run your example and tiler made the tiles without issue. preview.html looks good. This was with the update, maybe it will work for you now.

Session:

R version 3.5.0 (2018-04-23)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 16299)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    
attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     
other attached packages:
[1] tiler_0.2.9.9000
loaded via a namespace (and not attached):
[1] compiler_3.5.0

However, I also cannot make this display in Leaflet with the leaflet package.

After seeing your example, I decided to push some new convenience functions to the projections branch for pulling metadata out of tilemapresource.xml. You can update your leafletCRS call and everything should pass an identical check vs. the original written-out version.

E.g.,

meta <- tile_meta(tile_dir)
crs <- leafletCRS(crsClass = "L.Proj.CRS", code = "EPSG:3571",
                  proj4def = "+proj=laea +lat_0=90 +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
                  resolutions = meta$res, origin = meta$origin, bounds = meta$bounds)

[jmlondon]

Great! ... I'll try to update my example to have tiler handle all of the tiling. Maybe @tim-salabim or @schloerke could offer some guidance in getting things working in leaflet or mapview.

[jmlondon]

Ok ... so, everything appears to work (although, I still get a bunch of the "partial argument ..." warnings) and I no longer get any out of memory error.

But ... in the end, the tile directory is empty. So, it seems to do a bunch of work, but then not actually save the results to the tile directory. If I had to guess, the issue is likely some subtle difference in specification of paths within the system call on Windows vs macOS. I can try and do some additional troubleshooting along those lines, but it may not be until next week.

R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.4

[leonawicz]

I was hoping I'd be able to reproduce your warnings but no luck here and I don't have a Mac to test on. I am unsure where the along and along.with are coming from. It must be something more internal to what I wrote.

If I'm thinking about your tif correctly, no processing needs to be done so the vast majority of code here should be completely skipped. Your warnings appear to occur prior to the preparing for tiling statement.

Also nothing is jumping out at me that could be triggering those warnings earlier in the code lines for tile that run prior to calling .proj_check. Do you get those warnings under any other reproducible circumstances, maybe using a small tif, or single layer tif, WGS84 tif, anything else?

[jmlondon]

Those along and along.with warnings actually appear throughout the process. My initial googling hasn't produced any obvious hints.

I'll try some of the examples in the vignette and report back

[jmlondon]

Here's what I see from the first example in the vignette.

library(tiler)                                             
library(raster)                                            
#> Loading required package: sp
tile_dir <- file.path(tempdir(), "tiles")                  
map <- system.file("maps/map_wgs84.tif", package = "tiler")
(r <- raster(map))                                         
#> class       : RasterLayer 
#> dimensions  : 32, 71, 2272  (nrow, ncol, ncell)
#> resolution  : 0.8333333, 0.8333333  (x, y)
#> extent      : -125.0208, -65.85417, 23.27083, 49.9375  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#> data source : /Users/josh.london/Library/R/3.5/library/tiler/maps/map_wgs84.tif 
#> names       : map_wgs84 
#> values      : -0.7205295, 5.545086  (min, max)

tile(map, tile_dir, "0-3")                                 
#> Coloring raster...
#> Preparing for tiling...
#> Creating tiles. Please wait...
#> Creating tile viewer...
#> Complete.
## Warning messages:
## 1: In seq.default(along = band) :
##   partial argument match of 'along' to 'along.with'
## 2: In seq.default(along = band) :
##   partial argument match of 'along' to 'along.with'
## 3: In rep(strsplit(as.character(zoom), "-")[[1]], length = 2) :
##   partial argument match of 'length' to 'length.out'
list.files(tile_dir)                                       
#> character(0)

Session info

``` r devtools::session_info() #> ─ Session info ────────────────────────────────────────────────────────── #> setting value #> version R version 3.5.0 (2018-04-23) #> os macOS High Sierra 10.13.4 #> system x86_64, darwin15.6.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> tz America/Los_Angeles #> date 2018-06-15 #> #> ─ Packages ────────────────────────────────────────────────────────────── #> package * version date source #> assertthat 0.2.0 2017-04-11 CRAN (R 3.5.0) #> backports 1.1.2 2017-12-13 CRAN (R 3.5.0) #> callr 2.0.4 2018-05-15 CRAN (R 3.5.0) #> cli 1.0.0 2017-11-05 CRAN (R 3.5.0) #> clisymbols 1.2.0 2017-05-21 CRAN (R 3.5.0) #> crayon 1.3.4 2017-09-16 CRAN (R 3.5.0) #> debugme 1.1.0 2017-10-22 CRAN (R 3.5.0) #> desc 1.2.0 2018-05-01 cran (@1.2.0) #> devtools 1.13.5.9000 2018-05-08 Github (r-lib/devtools@13ee56b) #> digest 0.6.15 2018-01-28 CRAN (R 3.5.0) #> evaluate 0.10.1 2017-06-24 CRAN (R 3.5.0) #> htmltools 0.3.6 2017-04-28 CRAN (R 3.5.0) #> knitr 1.20 2018-02-20 CRAN (R 3.5.0) #> lattice 0.20-35 2017-03-25 CRAN (R 3.5.0) #> magrittr 1.5 2014-11-22 CRAN (R 3.5.0) #> memoise 1.1.0 2017-04-21 CRAN (R 3.5.0) #> pkgbuild 1.0.0 2018-05-08 Github (r-lib/pkgbuild@4ee12af) #> pkgload 1.0.0 2018-05-08 Github (r-lib/pkgload@35efedd) #> processx 3.1.0 2018-05-15 CRAN (R 3.5.0) #> R6 2.2.2 2017-06-17 CRAN (R 3.5.0) #> raster * 2.6-7 2017-11-13 CRAN (R 3.5.0) #> Rcpp 0.12.17 2018-05-18 cran (@0.12.17) #> rgdal 1.3-2 2018-06-08 CRAN (R 3.5.0) #> rlang 0.2.1 2018-05-30 CRAN (R 3.5.0) #> rmarkdown 1.9 2018-03-01 CRAN (R 3.5.0) #> rprojroot 1.3-2 2018-01-03 CRAN (R 3.5.0) #> sessioninfo 1.0.0 2017-06-21 cran (@1.0.0) #> sp * 1.3-1 2018-06-05 CRAN (R 3.5.0) #> stringi 1.2.2 2018-05-02 cran (@1.2.2) #> stringr 1.3.1 2018-05-10 cran (@1.3.1) #> testthat 2.0.0 2017-12-13 CRAN (R 3.5.0) #> tiler * 0.2.9.9000 2018-06-15 Github (leonawicz/tiler@942ea65) #> usethis 1.3.0 2018-02-24 CRAN (R 3.5.0) #> withr 2.1.2 2018-03-15 CRAN (R 3.5.0) #> yaml 2.1.19 2018-05-01 CRAN (R 3.5.0) ```

[jmlondon]

Ok, I think I've figured out what's wrong ... but, I don't understand what's going on

Bottom line, if I run the code from within R in a terminal session, it all works and the tiles are generated. But, if I run the code from within RStudio (both the release and the latest daily), the tiles are not generated. There must be some weird interaction between RStudio and the system call??

Bizarre

[jmlondon]

In further developments, if I open RStudio from the terminal via:

open -a RStudio

Then, the example code runs fine. There must be some environment variable or permissions setting that doesn't get picked up when RStudio is launched from the application icon. Annoying ... and, likely, beyond my skillset to come up with a fix.

[leonawicz]

Strange. The closest I've come to that on Windows is #5 I wonder if they are related in some way. I don't have answers to this either.

I have updated the branch, mostly a bunch of Windows system edge cases. It shouldn't affect anything for you. But do let me know if I manged to break anything in the process.

But aside from these Windows and Mac and RStudio things, it looks like we can continue focusing on the projections.

By the way, these tiles from your example took FOREVER to push to GitHub lol. Looks really good in a raw Leaflet widget. We just need to figure out these leaflet package specifications.

[leonawicz]

Closing since this issue is very old.