rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
543 stars 90 forks source link

Getting the wrong extention for a raster (bug) #1618

Closed bienflorencia closed 1 month ago

bienflorencia commented 1 month ago

I'm getting an error when trying to resample a raster layer. The error appears when I create the raster to use as a template. Here's a reprex:

library(geodata)
library(rnaturalearth)
library(sf)
sf_use_s2(FALSE) # switch spherical geometry off
library(terra)

peru <- ne_countries(country = "Peru", returnclass = "sf")
extent <- round(st_bbox(peru),0) + c(-1,-1,1,1) # Bounding box plus 1 degree

extent
# xmin ymin xmax ymax 
#  -82  -19  -68    1 

geodata_path('data')

tavg <- worldclim_country(var = "tavg", country = "PER")
tavg <- crop(tavg, extent)

r <- rast(res = 10000, # Target cell size
          ext = extent, # Extent of the raster
          crs = "epsg:4326") # Projection of our rasters

tavg <- resample(tavg, r, method = "bilinear")

tavg
# class       : SpatRaster 
# dimensions  : 1, 1, 12  (nrow, ncol, nlyr)
# resolution  : 10000, 10000  (x, y)
# extent      : -82, 9918, -19, 9981  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# source(s)   : memory
# names       : PER_w~avg_1, PER_w~avg_2, PER_w~avg_3, PER_w~avg_4, PER_w~avg_5, PER_w~avg_6, ... 
# min values  :         NaN,         NaN,         NaN,         NaN,         NaN,         NaN, ... 
# max values  :         NaN,         NaN,         NaN,         NaN,         NaN,         NaN, ... 

ext(r)
# SpatExtent : -82, 9918, -19, 9981 (xmin, xmax, ymin, ymax)

If I try with vect(peru) I have the same error.

ext(vect(peru))
# SpatExtent : -81.4109425523995, -68.6650797186896, -18.3479753557089, -0.0572054988648603 (xmin, xmax, ymin, ymax)

r <- rast(res = 10000, # Target cell size
          ext = ext(vect(peru)), # Extent of the raster
          crs = "epsg:4326") # Projection of our rasters

ext(r)
# SpatExtent : -81.4109425523995, 9918.5890574476, -18.3479753557089, 9981.65202464429 (xmin, xmax, ymin, ymax)

Thanks!

packageVersion("terra")
# [1] ‘1.7.78’
terra::gdal(lib="all")
#    gdal     proj     geos 
# "3.5.3"  "9.1.0" "3.11.0"
rhijmans commented 1 month ago

Your are doing the equivalent of

r <- rast(res=10000, ext=c(-82, -68, -19, 1), crs="+proj=longlat")

And that gives you

r
#class       : SpatRaster 
#dimensions  : 1, 1, 1  (nrow, ncol, nlyr)
#resolution  : 10000, 10000  (x, y)
#extent      : -82, 9918, -19, 9981  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 

That happens because we need to accommodate the extent you requested with the resolution. Your resolution is 1000 degrees! That is, of course, non-sensical. I assume that you want a resolution of approximately 1000 m (~ 30 seconds of a degree, or 1/120 degrees).

r <- rast(res=1/120, ext=c(-82, -68, -19, 1), crs="+proj=longlat")
#class       : SpatRaster 
#dimensions  : 2400, 1680, 1  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : -82, -68, -19, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
bienflorencia commented 1 month ago

Thank you, Robert! The resolution unit was intended for a layer with projected coordinates (in meters), and I mistakenly didn't change it to crs=4326 and latlong units.