Issue with projectRaster() when turning global lon/lat map into Mollweide projection

[RS-eco]

This code used to work fine in the past, but now throws out a very weird global map, when projecting a WGS84 global map to Mollweide projection. Any ideas what I am doing wrong here?

library(reprex)
library(rgeos)
#> Loading required package: sp
#> Warning in fun(libname, pkgname): rgeos: versions of GEOS runtime 3.10.2-CAPI-1.16.0
#> and GEOS at installation 3.8.0-CAPI-1.13.1differ
#> rgeos version: 0.5-9, (SVN revision 684)
#>  GEOS runtime version: 3.10.2-CAPI-1.16.0 
#>  Please note that rgeos will be retired by the end of 2023,
#> plan transition to sf functions using GEOS at your earliest convenience.
#>  GEOS using OverlayNG
#>  Linking to sp version: 1.4-7 
#>  Polygon checking: TRUE
library(rworldmap)
#> ### Welcome to rworldmap ###
#> For a short introduction type :   vignette('rworldmap')
library(raster)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

# Global grid
r_grid <- raster(nrows=360, ncols=720, xmn=-180, xmx=180, ymn=-90, ymx=90, 
                         crs="+proj=longlat + datum=WGS84", vals=1)

# Obtain world map
outline <- gPolygonize(gNode(as(getMap(resolution = "high"), 
                                              "SpatialLines")))
outline <- gUnaryUnion(outline)
outline <- st_as_sf(outline)

r_grid <- mask(r_grid, outline)
plot(r_grid)

r_grid_moll <- projectRaster(r_grid, crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain
#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain

#> Warning in x@ptr$project(y): GDAL Error 1: Point outside of projection domain
#> Warning in x@ptr$project(y): GDAL Error 1: Reprojection failed, err = 2050,
#> further errors will be suppressed on the transform object.
#> Warning in vals[i] <- value: number of items to replace is not a multiple of
#> replacement length
plot(r_grid_moll)

# Session info
sessionInfo()
#> R version 4.2.1 (2022-06-23)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] sf_1.0-8        raster_3.6-3    rworldmap_1.3-6 rgeos_0.5-9    
#> [5] sp_1.5-0        reprex_2.0.2   
#> 
#> loaded via a namespace (and not attached):
#>  [1] spam_2.9-1         tidyselect_1.2.0   terra_1.6-17       xfun_0.34         
#>  [5] lattice_0.20-45    colorspace_2.0-3   vctrs_0.5.0        generics_0.1.3    
#>  [9] htmltools_0.5.3    viridisLite_0.4.1  yaml_2.3.6         utf8_1.2.2        
#> [13] rlang_1.0.6        e1071_1.7-12       pillar_1.8.1       foreign_0.8-82    
#> [17] glue_1.6.2         withr_2.5.0        DBI_1.1.3          lifecycle_1.0.3   
#> [21] stringr_1.4.1      fields_14.1        dotCall64_1.0-2    munsell_0.5.0     
#> [25] gtable_0.3.1       codetools_0.2-18   evaluate_0.17      knitr_1.40        
#> [29] fastmap_1.1.0      maptools_1.1-5     class_7.3-20       fansi_1.0.3       
#> [33] highr_0.9          Rcpp_1.0.9         KernSmooth_2.23-20 classInt_0.4-8    
#> [37] scales_1.2.1       fs_1.5.2           gridExtra_2.3      ggplot2_3.3.6     
#> [41] digest_0.6.30      stringi_1.7.8      dplyr_1.0.10       grid_4.2.1        
#> [45] rgdal_1.5-32       cli_3.4.1          tools_4.2.1        magrittr_2.0.3    
#> [49] maps_3.4.0         proxy_0.4-27       tibble_3.1.8       rworldxtra_1.01   
#> [53] pkgconfig_2.0.3    assertthat_0.2.1   rmarkdown_2.17     rstudioapi_0.14   
#> [57] viridis_0.6.2      R6_2.5.1           units_0.8-0        compiler_4.2.1

[mdsumner]

I can't actually reproduce, I'm only on raster_3.5-29 and terra_1.6-7.

But generally use a target grid rather than a bare CRS - sometimes the hueristics to determine a grid are unhelpful. That might be a workaround for now.

ex <- c(-1, 1, -1/2, 1/2) * pi * 6378137
plot(projectRaster(r_grid, raster(extent(ex), res = 1e5, crs = "+proj=moll")))

[rhijmans]

Thank you for reporting this. This indeed fails with the CRAN version (@mdsumner is one version behind). The current version no longer uses "rgeos" and "rgdal"; the functionality these packages provided now comes from "terra". This transition introduced a few bugs.

This below (your example, simplified), now works for me with the development versions of both "raster" and "terra"

library(rworldmap)
library(raster)
world <- getMap(resolution = "low")
r <- rasterize(outline, raster())
m <- projectRaster(r, crs="+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
plot(m)