Is it possible to display a projected raster on a projected basemap in R/leaflet? The addRasterImage()
functions says that the projection needs to be in EPSG:3857. Can this not be changed by setting project = false
? I am able to display projected vector data on a projected basemap, but not raster ...
My attempt:
library(leaflet)
library(raster)
library(sf)
# Find location in northern Canada
ca_df <- data.frame(long = -114.3717401, lat = 62.4525548, name="Yellowknife", stringsAsFactors = F )
ca_pt <- st_as_sf(ca_df,coords = c("long", "lat"), crs = 4326)
# Project to Alaska Polar Stereographic
ca_pt_5936 <- as_Spatial(st_transform(ca_pt, 5936))@coords
# Create raster around point
r_5936 <- raster(
matrix(round(runif(100)), ncol = 10),
xmn = ca_pt_5936[[1]] - 50000, xmx = ca_pt_5936[[1]] + 50000,
ymn = ca_pt_5936[[2]] - 50000, ymx = ca_pt_5936[[2]] + 50000,
crs = "EPSG:5936"
)
# Project raster to Web Mercator (needed to get the extent in lat/long)
r_3857 <- projectRaster(r_5936, crs="EPSG:3857", method = "ngb")
# Prep for leaflet: https://github.com/rstudio/leaflet/issues/550
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)
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
)
# Map
leaflet(r_3857,
options= leafletOptions(
crs=epsg5936)) %>%
addTiles(urlTemplate = tile_url,
attribution = "Esri, DeLorme, GEBCO, NOAA NGDC, and other contributors",
options = tileOptions(minZoom = 0, maxZoom = 4)) %>%
addRasterImage(r_5936, project = F)
The output doesn't display the raster.