How to change SpatRast resolution using res and projectRaster in library(terra) correctly?-CodePudding

I have a raster file of sea surface temperatures from HadISST at 1x1 resolution. Previously I've re-grid to 0.5 x 0.5 using cdo(csl("remapbil","global_0.5"), inputfile, output file, debug=FALSE) from cdo via climateoperators but I'd like to explore how the workflow would be in terra using projectRaster. Example below using a HadISST nc file and rnaturalearth:

library(tidyverse)
library(ggplot2)
library(tidyterra)
library(terra)
library(tmaptools)
library(sf)
library(rnaturalearth)

### HadISST SST ndcf file (518kb)
if (!file.exists("HadISST_SST.nc.gz")) { 
     download.file("https://www.metoffice.gov.uk/hadobs/hadisst/data/HadISST1_SST_update.nc.gz","HadISST_SST.nc.gz")
     R.utils:::gunzip("HadISST_SST.nc.gz")
}

hadISST <- terra::rast('HadISST_SST.nc') 
hadISST_mean <- mean(hadISST) %>% clamp(lower=5, value=5)
hadISST_outlines <- mean(hadISST) %>% as.polygons(dissolve=FALSE)

### natural earth world coastline
worldrn_m <- ne_coastline(scale = "medium", returnclass = "sf") %>%  st_transform(crs = st_crs(hadISST))

ggplot()   theme_bw()   coord_fixed()   xlim(90,150)    ylim(-20,40)   
  geom_spatraster(data=hadISST_mean,  alpha=0.8)  
  geom_spatvector(data=hadISST_outlines, color="black", size=0.1, fill=NA,  alpha=0.5)   
  scale_fill_continuous(type = "viridis", limits=c(0, 35), breaks=seq(5,35,by=5))  
  geom_sf(data=worldrn_m)

projectRaster failed without a CRS using the following:

test <- projectRaster(hadISST_mean, res = 0.5, method="bilinear")

Error in rgdal::checkCRSArgs_ng(uprojargs = uprojargs, SRS_string = SRS_string,  : 
  is.character(SRS_string) is not TRUE

the crs I've used with sf previously throws the same error:

target_crs <- st_crs(" proj=longlat  x_0=0  y_0=0  lat_0=0  lon_0=180")
test <- projectRaster(hadISST_mean, res = 0.5, crs=target_crs, method="bilinear")

Error in rgdal::checkCRSArgs_ng(uprojargs = uprojargs, SRS_string = SRS_string,  : 
  is.character(SRS_string) is not TRUE

Can anyone advise where I'm going wrong here?

CodePudding user response：

There is no projectRaster method for a <SpatRaster>. You can use terra::project instead, either by providing a character representation of a coordinate reference system (see below), or (much better) by providing a raster template (see ?terra::project).

library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
pr <- project(r, "EPSG:2169", res=100)

And for a SpatVector

v <- vect(system.file("ex/lux.shp", package="terra"))
pv <- project(v, crs(pr))

plot(pr); lines(pv)

But if your goal is to change the resolution (as opposed to the coordinate reference system) project is not the appropriate method, and you should use disaggregate or resample instead. For example:

rr <- disagg(r, 2)

When asking questions, try to use a minimal, self-contained, example, as I show above; without loading unnecessary packages (never use library(tidyverse) to avoid asking us to install many packages, none of them relevant to the question at hand). Here is a reduced form of your example

library(terra)
library(geodata)
w <- geodata::world(path=".")
h <- rast(res=1) |> init("col") |> mask(w)
crs <- " proj=longlat  x_0=0  y_0=0  lat_0=0  lon_0=180"
test <- project(h, crs, res=0.5, method="bilinear")

I do not know why, but the projection is not working properly (nothing changes). The below does shift the location of the cropped area, but the values are gone.

hcrop <- crop(h, ext(90,150,-20,40)) 
test2 <- project(hcrop, crs, res=0.5, method="bilinear")
#class       : SpatRaster 
#dimensions  : 120, 120, 1  (nrow, ncol, nlyr)
#resolution  : 0.5, 0.5  (x, y)
#extent      : -90, -30, -20, 40  (xmin, xmax, ymin, ymax)
#coord. ref. :  proj=longlat  x_0=0  y_0=0  lat_0=0  lon_0=180 
#source      : memory 
#name        : lyr.1 
#min value   :   NaN 
#max value   :   NaN

The good news is that it seems that a better way to achieve what you are after goes like this:

x <- rotate(h, FALSE)

And you can change the resolution with:

y <- disagg(x, 2, method="bilinear")
plot(y)