|
| 1 | +#' --- |
| 2 | +#' title: "Polar Maps in Leaflet" |
| 3 | +#' author: "Bhaskar V. Karambelkar" |
| 4 | +#' --- |
| 5 | + |
| 6 | +library(leaflet) |
| 7 | + |
| 8 | +#' ## Artic Projections |
| 9 | + |
| 10 | +#' There is a [polarmap.js](http://webmap.arcticconnect.org/) |
| 11 | +#' leaflet plugin available, but that one is not easy to integrate in to the R package.<br/> |
| 12 | +#' But thankfully it does provide Tiles in different projections |
| 13 | +#' which can be used with Proj4Leaflet. |
| 14 | +#' In all it supports 6 projections and corresponding tile layers. |
| 15 | +#' <br/> |
| 16 | +#' The polarmap.js supports only Artic data, for Antartica see the end of this document. |
| 17 | + |
| 18 | + |
| 19 | +#' All these numbers and calculations come from the polarmap.js plugin, specifically from these files |
| 20 | +#' |
| 21 | +#' - http://webmap.arcticconnect.org/polarmap.js/dist/polarmap-src.js |
| 22 | +#' - http://webmap.arcticconnect.org/tiles.html |
| 23 | +#' - http://webmap.arcticconnect.org/usage.html |
| 24 | +#' |
| 25 | +extent <- 11000000 + 9036842.762 + 667 |
| 26 | +origin = c(-extent, extent) |
| 27 | +maxResolution <- ((extent - -extent) / 256) |
| 28 | +center <- c(90,0) |
| 29 | +defZoom <- 4 |
| 30 | +bounds <- list(c(-extent, extent),c(extent, -extent)) |
| 31 | +minZoom <- 0 |
| 32 | +maxZoom <- 18 |
| 33 | +resolutions <- purrr::map_dbl(minZoom:maxZoom,function(x) maxResolution/(2^x)) |
| 34 | + |
| 35 | +# 6 Projection EPSG Codes |
| 36 | +projections <- c('3571', '3572', '3573', '3574', '3575', '3576') |
| 37 | +# Corresponding proj4defs codes for each projection |
| 38 | +proj4defs <- list( |
| 39 | + '3571' = '+proj=laea +lat_0=90 +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', |
| 40 | + '3572' = '+proj=laea +lat_0=90 +lon_0=-150 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', |
| 41 | + '3573' = '+proj=laea +lat_0=90 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', |
| 42 | + '3574' = '+proj=laea +lat_0=90 +lon_0=-40 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', |
| 43 | + '3575' = '+proj=laea +lat_0=90 +lon_0=10 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', |
| 44 | + '3576' = '+proj=laea +lat_0=90 +lon_0=90 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs' |
| 45 | +) |
| 46 | + |
| 47 | +# create a CRS instance for each projection |
| 48 | +crses <- purrr::map(projections, function(code) { |
| 49 | + leafletCRS( |
| 50 | + crsClass = 'L.Proj.CRS', |
| 51 | + code = sprintf("EPSG:%s",code), |
| 52 | + proj4def = proj4defs[[code]], |
| 53 | + origin = origin, |
| 54 | + resolutions = resolutions, |
| 55 | + bounds = bounds |
| 56 | + ) |
| 57 | +}) |
| 58 | + |
| 59 | +# Tile URL Template for each projection |
| 60 | +tileURLtemplates <- purrr::map(projections, function(code) { |
| 61 | + sprintf('http://{s}.tiles.arcticconnect.org/osm_%s/{z}/{x}/{y}.png', |
| 62 | + code) |
| 63 | +}) |
| 64 | + |
| 65 | +# We can't add all 6 tiles to our leaflet map, |
| 66 | +# because each one is in a different projection, |
| 67 | +# and you can have only one projection per map in Leaflet. |
| 68 | +# So we create 6 maps. |
| 69 | +polarmaps <- purrr::map2(crses, tileURLtemplates, |
| 70 | + function(crs, tileURLTemplate) { |
| 71 | + leaflet(options= list(crs=crs, zoom=defZoom, center = center)) %>% |
| 72 | + addTiles(urlTemplate = tileURLTemplate, |
| 73 | + attribution = "Map © ArcticConnect. Data © OpenStreetMap contributors", |
| 74 | + options = tileOptions(subdomains = "abc", noWrap = TRUE, |
| 75 | + continuousWorld = FALSE)) |
| 76 | + }) |
| 77 | + |
| 78 | +#' #### EPSG:3571 |
| 79 | +polarmaps[[1]] %>% |
| 80 | + addGraticule() |
| 81 | + |
| 82 | +#' #### EPSG:3572 |
| 83 | +polarmaps[[2]] |
| 84 | + |
| 85 | +#' #### EPSG:3573 |
| 86 | +polarmaps[[3]] |
| 87 | + |
| 88 | +#' #### EPSG:3574 |
| 89 | +polarmaps[[4]] |
| 90 | + |
| 91 | +#' #### EPSG:3575 |
| 92 | +polarmaps[[5]] |
| 93 | + |
| 94 | +#' #### EPSG:3576 |
| 95 | +polarmaps[[6]] |
| 96 | + |
| 97 | +#' ## Antartica |
| 98 | +#' ### NOT YET WORKING |
| 99 | +#' Code adapted from |
| 100 | +#' https://github.com/nasa-gibs/gibs-web-examples/blob/release/examples/leaflet/antarctic-epsg3031.js <br/> |
| 101 | +#' This should work in theory but not working yet. |
| 102 | + |
| 103 | +resolutions <- c(8192, 4096, 2048, 1024, 512, 256) |
| 104 | +center <- c(-90,0) |
| 105 | +zoom <- 0 |
| 106 | +maxZoom <- 5 |
| 107 | + |
| 108 | +crsAntartica <- leafletCRS( |
| 109 | + crsClass = 'L.Proj.CRS', |
| 110 | + code = 'EPSG:3031', |
| 111 | + proj4def = '+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs', |
| 112 | + resolutions = resolutions, |
| 113 | + origin = c(-4194304, 4194304), |
| 114 | + bounds = list( c(-4194304, -4194304), c(4194304, 4194304) ) |
| 115 | +) |
| 116 | + |
| 117 | +antarticaTilesURL <- "//map1{s}.vis.earthdata.nasa.gov/wmts-antarctic/MODIS_Aqua_CorrectedReflectance_TrueColor/default/2014-12-01/EPSG3031_250m/{z}/{y}/{x}.jpg" |
| 118 | + |
| 119 | +leaflet(options= list( |
| 120 | + crs=crsAntartica, |
| 121 | + zoom=zoom, maxZoom=maxZoom, center= center, |
| 122 | + continuousWorld = TRUE, worldCopyJump = FALSE)) %>% |
| 123 | + setView(-90,0,0) %>% |
| 124 | + addTiles(urlTemplate = antarticaTilesURL, |
| 125 | + layerId = "antartica_tiles", |
| 126 | + attribution = "<a href='https://earthdata.nasa.gov/gibs'> NASA EOSDIS GIBS</a> <a href='https://github.com/nasa-gibs/web-examples/blob/release/leaflet/js/antarctic-epsg3031.js'> View Source </a>", |
| 127 | + options = tileOptions( |
| 128 | + tileSize =512, |
| 129 | + subdomains = "abc", |
| 130 | + noWrap = TRUE, |
| 131 | + continuousWorld = TRUE, |
| 132 | + format = "image%2Fjpeg" |
| 133 | + )) %>% |
| 134 | + htmlwidgets::onRender( |
| 135 | + "function(el, t){ |
| 136 | + var myMap = this; |
| 137 | + debugger; |
| 138 | + var tileLayer = myMap.layerManager._byLayerId['tile\\nantartica_tiles']; |
| 139 | +
|
| 140 | + // HACK: BEGIN |
| 141 | + // Leaflet does not yet handle these kind of projections nicely. Monkey |
| 142 | + // patch the getTileUrl function to ensure requests are within |
| 143 | + // tile matrix set boundaries. |
| 144 | + var superGetTileUrl = tileLayer.getTileUrl; |
| 145 | +
|
| 146 | + tileLayer.getTileUrl = function(coords) { |
| 147 | + debugger; |
| 148 | + var max = Math.pow(2, tileLayer._getZoomForUrl() + 1); |
| 149 | + if ( coords.x < 0 ) { return ''; } |
| 150 | + if ( coords.y < 0 ) { return ''; } |
| 151 | + if ( coords.x >= max ) { return ''; } |
| 152 | + if ( coords.y >= max ) { return ''; } |
| 153 | + return superGetTileUrl.call(tileLayer, coords); |
| 154 | + }; |
| 155 | + // HACK: END |
| 156 | + }") |
0 commit comments