2

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.

enter image description here

0

Browse other questions tagged or ask your own question.