Shift in Raster image in R Leaflet

Question

I am trying to plot a raster image on a leaflet map - the raster I am using is downloaded from Environment Agency and the CRS of the raster is British National Grid (EPSG:27700). You can find the file here

rast <- raster('test.tif')
proj4string(rast) #confirming crs
[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_d

Plotting it using leaflet

leaflet()%>% 
  addTiles() %>% 
  addRasterImage(rast)

I don't understand why there is a shift here (you can see it from the building footprint on the map and the raster). In the documentation for addRasterImage, it says

By default, the addRasterImage function will project the RasterLayer x to EPSG:3857 using the raster package's projectRaster function.

So it should basically work with any CRS. When I plot rasters with a crs EPSG:4326 there are no shifts and everything looks fine.

Answer 1

I think this is related to problems with PROJ4 strings, and hence the move to WKT projection strings in modern PROJ systems.

The quick fix is to reset the projection to epsg:27700 via a CRS call:

> rast <- raster('test.tif')
> proj4string(rast) = CRS("+init=epsg:27700")

mapview(rast) now puts it in place on a web map.

I'm not sure what the difference is between CRS("+init=epsg:27700"), which produces a warning:

CRS arguments:
 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
+y_0=-100000 +ellps=airy +units=m +no_defs 
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum OSGB_1936 in Proj4 definition

...and the proj4 string from the raster:

> proj4string( raster('test.tif') )
[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"

Which produced a different warning when I tried to work with it:

Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4 definition

Here's the code I'm using and the maps:

library(raster)
rast <- raster('test.tif')
library(mapview)
mapview(rast)

proj4string(rast) = CRS("+init=epsg:27700")
mapview(rast)

I think its because the proj4 string doesn't well-define a datum, but the epsg code 27700 does. Using sf package I get:

> st_crs(proj4string( raster('test.tif') ))
Coordinate Reference System:
  User input: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["Unknown based on Airy 1830 ellipsoid",
            ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],

which has an unknown datum, but

> st_crs(27700)
Coordinate Reference System:
  User input: EPSG:27700 
  wkt:
PROJCRS["OSGB 1936 / British National Grid",
    BASEGEOGCRS["OSGB 1936",
        DATUM["OSGB 1936",
            ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,

has a well-known OSGB36 datum. This has a shift you can see with proj on the command line:

     OSGB36 airy         towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894

and I think that towgs84 is the practical difference here.

QGIS loads the tiff in the correct place with no need for adjustment, so it may have a better idea of how to shift these things correctly...