rstudio / leaflet

R Interface to Leaflet Maps
http://rstudio.github.io/leaflet/
Other
805 stars 509 forks source link

support for ESRI Arctic Ocean Base map (polar stereo projection) #550

Open jmlondon opened 6 years ago

jmlondon commented 6 years ago

I have been searching for a viable solution that would allow use of leaflet (and, also, mapview) with points, lines, polygons and rasters that cross the 180 longitude line. I have long hoped for a reliable, simple solution based on web mercator. But, I've yet to find anything satisfactory.

Leaflet's support of additional projections is a great solution except for the lack of tile options available. That said, ESRI has recently released Arctic Ocean Basemap which is in the Alaska Polar Stereographic projection.

I've managed to sort out all the parameter details for the tile_url and leafletCRS function. Here, I present an example that others who are interested in the US Arctic may find useful.

I would also be interested in how I might contribute this to either leaflet or, maybe leaflet-providers or leaflet-extras, so the full configuration doesn't have to be specified each time. I imagine something similar could be specified for ArcticConnect and other tile sources (e.g. Antarctica).

library(leaflet)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
library(sf)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
#> Linking to GEOS 3.6.2, GDAL 2.2.4, proj.4 5.0.1

track_line <-                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
'LINESTRING (-169.4444 62.57549, -169.8999 62.31869, -170.4738 62.16436, -170.6529 61.4344, -171.7103 60.5111, -173.5295 60.78438, -175.2325 61.55942, -176.0107 61.88525, -177.1819 61.43621, -177.3077 61.25797, -177.3493 61.98553, -176.4385 62.26207, -177.3228 62.73083, -179.2457 62.84639, 179.7168 62.67369, 179.8549 62.58532, 179.8978 62.60145, -179.7893 62.64932, 179.9897 62.37538, -179.6393 62.66016, -179.8876 62.60982, -179.68 62.75012, 179.8264 62.78815, -179.789 62.72092, 179.6381 62.48176)'

track_line <- sf::st_as_sfc(track_line) %>%                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
sf::st_set_crs(4326)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  

tile_url <- 'https://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer/tile/{z}/{y}/{x}.png'                                                                                                                                                                                                                                                                                                                                                                                           
origin <- c(-2.8567784109255e+07, 3.2567784109255e+07)                                                                                                                                                                                                                                                                                                                                                                                                                                                                
minZoom <- 0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
maxZoom <- 24                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
resolutions <- c(                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
238810.813354,119405.406677, 59702.7033384999, 29851.3516692501,14925.675834625,
7462.83791731252,3731.41895865639, 1865.70947932806,932.854739664032,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
466.427369832148, 233.213684916074, 116.60684245803701, 58.30342122888621,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
29.151710614575396, 14.5758553072877, 7.28792765351156, 3.64396382688807,
1.82198191331174, 0.910990956788164, 0.45549547826179, 0.227747739130895,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
0.113873869697739, 0.05693693484887, 0.028468467424435                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     

epsg5936 <- leafletCRS(                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
crsClass = 'L.Proj.CRS',                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              
code = 'EPSG:5936',                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
proj4def = '+proj=stere +lat_0=90 +lat_ts=90 +lon_0=-150 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs',                                                                                                                                                                                                                                                                                                                                                                                          
origin = origin,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
resolutions = resolutions                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     

leaflet(track_line,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
options= leafletOptions(                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              
crs=epsg5936)) %>%                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
addTiles(urlTemplate = tile_url,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
attribution = "Esri, DeLorme, GEBCO, NOAA NGDC, and other contributors",                                                                                                                                                                                                                                                                                                                                                                                                                                              
options = tileOptions(minZoom = 0, maxZoom = 24)) %>%                                                                                                                                                                                                                                                                                                                                                                                                                                                                 
addPolylines()                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        
schloerke commented 6 years ago
... cross the 180 longitude line. I have long hoped for a reliable, simple solution based
on web mercator. But, I've yet to find anything satisfactory.

I believe this is still the case for standard leaflet. I found a plugin here: https://github.com/briannaAndCo/Leaflet.Antimeridian . This would need a PR to this repo to be incorporated.

ESRI has recently released Arctic Ocean Basemap which is in the Alaska Polar 
Stereographic projection.  I've managed to sort out all the parameter details for the 
`tile_url` and `leafletCRS()` function

There is a third package: leaflet.esri. It contains a method called addEsriTiledMapLayer(). You can use it like:

leaflet() %>% 
  leaflet.esri::addEsriTiledMapLayer(
    "https://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer"
  )

The projection is still off, so you will still have to use the leafletCRS() method as you did.

I would also be interested in how I might contribute this

leaflet-providers has ESRI maps and would be a good place to submit a PR to add the map location. Eventually, it will trickle into leaflet R package.

I believe a leaflet.providers package needs to be made due to how quickly providers are added.