Support Rendering Raster in Original Projection

[DanielJDufour]

Is your feature request related to a problem? Please describe. People want to see the raster in the original projection on a Leaflet map.

Describe the solution you'd like People want to be able to load the raster and also change the map projection to that of the raster.

Describe alternatives you've considered

• need to look at Leaflet plugins for this

Additional context https://twitter.com/TimSalabim3/status/1319967430064693248

[mdsumner]

nt_20201024_f18_nrt_s.zip

That is NSIDC 25km daily sea ice concentration, converted from the raw binary file to GeoTIFF. in R it looks like

class      : RasterStack 
dimensions : 332, 316, 104912, 1  (nrow, ncol, ncell, nlayers)
resolution : 25000, 25000  (x, y)
extent     : -3950000, 3950000, -3950000, 4350000  (xmin, xmax, ymin, ymax)
crs        : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs 
names      : nt_20201024_f18_nrt_s.bin 
min values :                         0 
max values :                       100 
time        : 2020-10-24 

source file from ftp server at

sidads.colorado.edu/pub/DATASETS/nsidc0081_nrt_nasateam_seaice/south/nt_20201024_f18_nrt_s.bin

NSIDC SMMR-SSM/I Nasateam sea ice concentration

Passive microwave estimates of sea ice concentration at 25km spatial resolution. Daily and monthly resolution, available from 1-Oct-1978 to present. Data undergo a quality checking process and are updated annually.

http://nsidc.org/data/nsidc-0051.html

[mdsumner]

Can easily do any kind of example, bathymetry topography etc. (I'm focussed on something else rn otherwise I'd look more closely) HTH Here's a bigger one, 8m topography extract of the Vestfold Hills from REMA 8m - aggregated by mean to 32m

vestfold.zip

Reference Elevation Model of Antarctica mosaic tiles
The Reference Elevation Model of Antarctica (REMA) is a high resolution, time-stamped digital surface model of Antarctica at 8-meter spatial resolution. REMA is constructed from hundreds of thousands of individual stereoscopic Digital Elevation Models (DEM) extracted from pairs of submeter (0.32 to 0.5 m) resolution DigitalGlobe satellite imagery. Version 1 of REMA includes approximately 98% of the contiguous continental landmass extending to maximum of roughly 88 degrees S. Output DEM raster files are being made available as both 'strip' files as they are output directly from SETSM that preserve the original source material temporal resolution, as well as mosaic tiles that are compiled from multiple strips that have been co-registered, blended, and feathered to reduce edge-matching artifacts.

class      : RasterLayer 
dimensions : 880, 950, 836000  (nrow, ncol, ncell)
resolution : 32, 32  (x, y)
extent     : 2301440, 2331840, 467528, 495688  (xmin, xmax, ymin, ymax)
crs        : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source     : vestfold.tif 
names      : vestfold 
values     : 18.10381, 332.6073  (min, max)

https://www.pgc.umn.edu/data/rema/

[tim-salabim]

Thanks Mike!

[DanielJDufour]

Thank you @mdsumner . The examples were very helpful (and really quite fascinating).

@tim-salabim , I published a new version of georaster-layer-for-leaflet that works with L.CRS.Simple. It displays the given raster without any reprojection and still uses nearest neighbor sampling. Because the map is technically using the "simple" projection and not the projection of the original raster, you won't be able to easily add other layers on top of the raster.

Future work (which would be a much bigger lift and won't happen for a while) includes:

• supporting maps in arbitrary projections
• adding different sampling algorithms - bilinear, etc.

Instructions: You don't have to do much to display a raster with the simple projection. There's actually no new parameters you have to use or change for creating a GeoRasterLayer. You just have to create the Leaflet map with the simple CRS like:

var map = L.map('map', {
    crs: L.CRS.Simple
});

And, of course, make sure to run

map.fitBounds(layer.getBounds());

You can view the source code for a test of the sea-ice around Antartica here:

• https://github.com/GeoTIFF/georaster-layer-for-leaflet/blob/master/tests/simple.html

I haven't update the official documentation for this library because I wanted to first see how it works for you. Would you be able to give it a try and let me know how it goes?

Thank you!!

[mdsumner]

coooool! Glad it's helpful, very happy to make any other examples in any projection (I can do imagery as well there are some interesting maps as well as sat imagery etc).

[tim-salabim]

Thanks @DanielJDufour ! This is fantastic!

Because the map is technically using the "simple" projection and not the projection of the original raster, you won't be able to easily add other layers on top of the raster.

I may be mistaken, but as long as the additional layers have identical projection to the raster, overlays should work. I will certainly try as easy adding of layers is one of the key objectives of my efforts in this realm.

Future work (which would be a much bigger lift and won't happen for a while) includes:

• supporting maps in arbitrary projections
• adding different sampling algorithms - bilinear, etc.

Ambitious goals indeed. If I might ask for one more, which will hopefully be much easier to implement:

• have an option to turn off resampling.

There are different reasons for this:

• distortion of the raster values when you really want to check the true cell values - https://twitter.com/MichaelDorman84/status/1320821770933456897
• for small (low resolution) rasters it can lead to unexpected results (and doesn't really provide any performance benefit I assume) - https://github.com/r-spatial/mapview/issues/339

I will try to get the new version into leafem soon (probably over the weekend).

Again, thanks for a great plugin!!

[tim-salabim]

@DanielJDufour I just updated to the latest georaster-layer-for-leaflet version and can confirm that it works in principle with L.CRS.Simple. However, there are two issues:

1. I cannot zoom out beyond a zoom of 0, i.e. at zoom -1 the image disappears.
2. When overlaying vector data (e.g. markers) in the same projection within the extent of the raster, they appear somewhere else.

From my experience, it seems that the vector data are indeed adhering to the projection extent (whatever that will be in the native, non-webmercator projection), whereas the raster gets relocated so that the lower left corner is at 0/0. Is that correct?

Is there any way that we could have the raster display in L.CRS.Simple behave the same way the vector data is rendered? I.e. within their natural bounds as defined by their extent?

Here's the proof that it works:

[DanielJDufour]

Hi, @tim-salabim . Thank you for the feedback. It's very helpful. I'll try to address these issues.

Would you be able to share the code you wrote to test When overlaying vector data (e.g. markers) in the same projection within the extent of the raster, they appear somewhere else., so I can explore the issue more thoroughly.

I'll also work on adding in support for turning off sampling.

[tim-salabim]

lcrs_simple_overlay.zip

@DanielJDufour do you know R? All this is in R via a package called htmlwidgets which acts like a translator between R objects and JS to produce valid hmtl. The R code to produce the attached map is

library(leaflet)
library(leafem)
library(stars)
library(sf)

# tif = system.file("tif/L7_ETMs.tif", package = "stars")
tif = "/home/timpanse/Downloads/nt_20201024_f18_nrt_s.tif"
x1 = read_stars(tif)

pts = st_sample(st_as_sfc(st_bbox(x1)), 100)

pal = hcl.colors(256, "viridis")
# ndvi = function(rst) (rst[4] - rst[3]) / (rst[4] + rst[3])

m = leaflet(
  options = leaflet::leafletOptions(
    minZoom = -1000
    , maxZoom = 52
    , crs = leafletCRS(crsClass = "L.CRS.Simple")
  )
) %>%
  # addProviderTiles("Esri.WorldImagery") %>%
  addGeotiff(
    file = tif
    # , arith = ndvi
    , project = FALSE
    # , bands = c(3)
    , resolution = 500
    , opacity = 1
    , colorOptions = colorOptions(
      palette = pal
      , na.color = "magenta"
      , domain = c(0, 100)
    )
    , pixelValuesToColorFn = NULL #myCustomJSFunc
    , group = "raster"
  ) %>%
  addCircleMarkers(data = pts, group = "points") %>%
  addLayersControl(overlayGroups = "raster") %>%
  addMouseCoordinates() %>%
  addHomeButton(raster::extent(c(-1, 1000, -1, 1000)), group = "raster") %>%
  addHomeButton(st_bbox(pts), group = "points")

m

mapview::mapshot(m, url = "/home/timpanse/Desktop/lcrs_simple_overlay.html")

I hope that the attached map is useful. It should all be there.

I think the crux of the issue is whether the raster overlay in L.CRS.Simple has to be in pixel dimensions (0-X/0-X) or if there is a way to display it in it's native coordinates (like the points in the attached example).