Search code examples
rrastertmap

Raster is drawn only partly in tmap


Some of my maps that contain rasters get drawn partly only. The southern part of the map is cut off. It looks like some memory limit, but I have no idea if it is so. Another option is something concerning reprojection and bbox... Nevertheless, other rasters of the same crs work well in this project. Final map cut off in the south

library(tmap)
library(tmaptools)
library(sf)
library(raster)
library(maps)

raster.file <-   tempfile()
download.file("https://github.com/andliszmmu/storage/raw/main/copernicus_2000.tif", raster.file, mode="wb")

ocean.pal   <- '#afdff9'
prj         <- CRS("+proj=eqdc +lat_0=57 +lon_0=46.8 +lat_1=69.4367 +lat_2=51.3043 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs")
map_bbox    <- st_bbox(c(xmin = -1700000, ymin = -1780000, xmax = 1040000, ymax = 2330000), crs = prj)

land        <- st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))
land        <- land[land$ID=='Russia',]
land        <- st_transform(land, prj)

model        <- raster(paste0(raster.file))

  tm =  
    tm_shape(land, bbox = map_bbox) +
    tm_fill(col = ocean.pal) +
    tm_shape(model)+
    tm_raster(style='fixed',
              breaks = c(0, 35,45,55,67,75,85,95,112, 200),
              palette=c('lightgreen', 'yellow', 'red', 'grey', 'white', 'blue', '#f21deb', 'darkgreen', 'green'),
              labels=c("a","b","c","d","e","f","g","h","i"),
              alpha = 1,
              #n=5,
              stretch.palette = FALSE,
              title = 'Landscape classification',
              legend.show = TRUE,
              legend.reverse = TRUE,
              legend.z=0) +
    tm_layout(frame = T,
              frame.lwd = 2,
              inner.margins = c(0,0,0,0),
              outer.margins = c(0.07, 0.127, 0.07, 0.112),
              asp = 39/56,
              legend.frame = T,
              legend.format = list(text.separator = "–"),
              legend.frame.lwd = 1,
              legend.outside = F,
              legend.position = c(0.001, 0.999),
              legend.just = c("left", "top"),
              legend.width = -0.32,
              legend.height = 0.9,
              legend.text.size = 0.7,
              legend.title.size = 1.1,
              legend.bg.color = "white",
              title.size = 0.7,
              title.bg.color = 'white',
              title.position = c('left', 'bottom')) +
    tmap_options(check.and.fix = TRUE,
                 max.raster = c(plot = 1e8, view = 1e8))
  
tm
    

I checked documentation and the web, but had not found solution.
Parameters of this raster are:
class : RasterLayer
dimensions : 2000, 1441, 2882000 (nrow, ncol, ncell)
resolution : 0.0343, 0.0179 (x, y)
extent : 19.64043, 69.06673, 41.19558, 76.99558 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : copernicus_2000.tif
names : copernicus_2000
values : 20, 200 (min, max)
Raster downsampling (resolution <- 0.0343, 0.0343) did not help.


Solution

  • Those non-conformal projections can be difficult to negotiate. Curious that none of your other data were clipped but I'm assuming none extend as far south as the problem raster?

    Either way, here is a solution. It involves raster::extend()ing the southern extent of your raster:

    library(tmap)
    library(tmaptools)
    library(sf)
    library(raster)
    library(maps)
    
    # PROJ4 string
    prj <- "+proj=eqdc +lat_0=57 +lon_0=46.8 +lat_1=69.4367 +lat_2=51.3043 +x_0=0 +y_0=0 +ellps=sphere +units=m +no_defs"
    
    # Load .tif as raster from working directory, previously downloaded from
    # https://github.com/andliszmmu/storage/raw/main/copernicus_2000.tif
    r <- raster("data/copernicus_2000.tif")
    
    # Return extent of raster
    st_bbox(r)
    #     xmin     ymin     xmax     ymax 
    # 19.64043 41.19558 69.06673 76.99558 
    
    # Extend the southern extent of raster e.g. change 41.19558 to 30
    r <- extend(r, extent(19.64043, 69.06673, 30, 76.99558)) # xmin xmax ymin ymax 
    
    map_bbox <- st_bbox(c(xmin = -1700000, ymin = -1780000, xmax = 1040000, ymax = 2330000), crs = prj)
    
    ocean.pal <- '#afdff9'
    land <- st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))
    land <- land[land$ID=='Russia',]
    land <- st_transform(land, prj)
    
    tm =
      tm_shape(land, bbox = map_bbox) +
      tm_fill(col = ocean.pal) +
      tm_shape(r)+
      tm_raster(style = 'fixed',
                breaks = c(0, 35,45,55,67,75,85,95,112, 200),
                palette=c('lightgreen', 'yellow', 'red', 'grey', 'white', 'blue', '#f21deb', 'darkgreen', 'green'),
                labels=c("a","b","c","d","e","f","g","h","i"),
                alpha = 1,
                #n=5,
                stretch.palette = FALSE,
                title = 'Landscape classification',
                legend.show = TRUE,
                legend.reverse = TRUE,
                legend.z=0) +
      tm_layout(frame = T,
                frame.lwd = 2,
                inner.margins = c(0,0,0,0),
                outer.margins = c(0.07, 0.127, 0.07, 0.112),
                asp = 39/56,
                legend.frame = T,
                legend.format = list(text.separator = "–"),
                legend.frame.lwd = 1,
                legend.outside = F,
                legend.position = c(0.001, 0.999),
                legend.just = c("left", "top"),
                legend.width = -0.32,
                legend.height = 0.9,
                legend.text.size = 0.7,
                legend.title.size = 1.1,
                legend.bg.color = "white",
                title.size = 0.7,
                title.bg.color = 'white',
                title.position = c('left', 'bottom')) +
      tmap_options(check.and.fix = TRUE,
                   max.raster = c(plot = 1e8, view = 1e8))
    
    tm
    

    1