rggplot2r-rastermap-projections

R ggplot plotting map raster with rounded shape - How to remove data outside projected area?


I am trying to plot a raster in a projected in a coordinated system which follows the curvature of the earth like most projections that are not WGS84. The problem is that the places were the globe wraps around the data should not be plotted outside the globe. I realize that ggplot cannot do a rounded/elliptical plot but how do I mask or remove automatically the data outside the globe? I have to plot more than 100 maps and I can't do this manually especially if I want to change to a different projection.

There's an answer here but it's hackish and doesn't seem to apply to every case, is there function or package that deals with this problem? I don't think R users only plot maps in WGS84? I am attaching a file and code to quickly plot the map. I cannot use xlim because it would cut some parts of the map since the borders are not straight.

#netcdf file
https://ufile.io/fy08x33d
library(terra);library(tidyterra)
r=rast('Beck_KG_V1_present_0p5.tif')
#background map
r[r==0]=NA
ggplot() +geom_spatraster(data=r)+scale_fill_viridis_c(na.value='transparent') +coord_sf(crs=st_crs("+proj=hatano"),expand=FALSE)

enter image description here


Solution

  • With these data

    library(terra)
    library(tidyterra)
    r1 <- rast('Beck_KG_V1_present_0p5.tif')
    r <- subst(r1, 0, NA)
    

    You can do

    library(ggplot2)
    p <- project(r, method="near", "+proj=hatano", mask=TRUE)
    ggplot() +geom_spatraster(data=p)+scale_fill_viridis_c(na.value='transparent') 
    

    enter image description here

    And here are two alternatives with base plot

    First with your own color palette and a legend

    library(viridis)
    g <- graticule(60, 45, "+proj=hatano")
    plot(g, background="azure", mar=c(.2,.2,.2,4), lab.cex=0.5, col="light gray")
    plot(p, add=TRUE, axes=FALSE, plg=list(shrink=.8), col=viridis(25))
    

    enter image description here

    With the colors that came with the file:

    coltab(p) <- coltab(r1)
    plot(g, background="azure", mar=.5, lab.cex=0.5, col="light gray")
    plot(p, add=TRUE, axes=FALSE, col=viridis(25))
    

    enter image description here