CODEWITHSUNDEEP
rrasterresolution
perep1972
perep1972

Reputation: 147

R language, unreasonable (?) raster resolutions

Sorry for the very stupid question, but I'm really stuck here... I need to create a Digital Elevation Model for my study area. For this, I downloaded an SRTM (1 arc-seg resolution, freely available from the net) image, which comprises a region wider than my area of interest. The original raster has these characteristics:

class      : RasterLayer 
dimensions : 3601, 3601, 12967201  (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778  (x, y)
extent     : -45.00014, -43.99986, -22.00014, -20.99986  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : s22_w045_1arc_v3.tif 
names      : s22_w045_1arc_v3 
values     : -32768, 32767  (min, max)

I need to (1) increase the resolution (initially of 30.75662 * 28.68392 m) to 1 * 1 m (that is, I really do not care about the exactitude of the elevations) and (2) crop a squared area of 2000 * 2000 m centered at a given coordinate. So, the first step I'm following is to re-project to UTM:

projection(r) <- "+proj=utm +zone=23 +datum=WGS84"

But the resolution units do not change after that:

class      : RasterLayer 
dimensions : 3601, 3601, 12967201  (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778  (x, y)
extent     : -45.00014, -43.99986, -22.00014, -20.99986  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs 
source     : s22_w045_1arc_v3.tif 
names      : s22_w045_1arc_v3 
values     : -32768, 32767  (min, max)

If I try to set the resolutions in meters manually, then generates an empty raster. Can anybody be so kind as to throw some light on me here?

Upvotes: 1

Views: 78

Answers (1)

dieghernan
dieghernan

Reputation: 3402

You are changing (i.e. overwriting) the CRS, not projecting the raster. Usually, it is recommended to create a template raster with the CRS and the resolution you need and reproject the raster using this template.

See here an example, I am switching to terra for the analysis since it is a newer and faster package, but I would show also how to convert it back to raster format:

library(raster)
#> Loading required package: sp

# Faking your data
r <- raster(
  nrows = 3601, ncols = 3601,
  ext = extent(c(-45.00014, -43.99986, -22.00014, -20.99986))
)



values(r) <- seq(-32768, 32767, length.out = ncell(r))
r
#> class      : RasterLayer 
#> dimensions : 3601, 3601, 12967201  (nrow, ncol, ncell)
#> resolution : 0.0002777784, 0.0002777784  (x, y)
#> extent     : -45.00014, -43.99986, -22.00014, -20.99986  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : memory
#> names      : layer 
#> values     : -32768, 32767  (min, max)
plot(r)


# End of faking data

# Change to terra, much faster
library(terra)
#> terra 1.5.21
r_terra <- terra::rast(r)
template <- terra::project(r_terra, "+proj=utm +zone=23 +datum=WGS84")

# Change to the desired res
res(template) <- c(2000, 2000)


# Reproject
r_terra_reproj <- terra::project(r_terra, template)

r_terra_reproj
#> class       : SpatRaster 
#> dimensions  : 56, 52, 1  (nrow, ncol, nlyr)
#> resolution  : 2000, 2000  (x, y)
#> extent      : 499985.4, 603985.4, -2433195, -2321195  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs 
#> source      : memory 
#> name        :     layer 
#> min value   : -32371.95 
#> max value   :  32182.81


terra::plot(r_terra_reproj)



# Back to RasterLayer

r_reproj <- raster(r_terra_reproj)

r_reproj
#> class      : RasterLayer 
#> dimensions : 56, 52, 2912  (nrow, ncol, ncell)
#> resolution : 2000, 2000  (x, y)
#> extent     : 499985.4, 603985.4, -2433195, -2321195  (xmin, xmax, ymin, ymax)
#> crs        : +proj=utm +zone=23 +datum=WGS84 +units=m +no_defs 
#> source     : memory
#> names      : layer 
#> values     : -32371.95, 32182.81  (min, max)

Created on 2022-06-10 by the reprex package (v2.0.1)

Upvotes: 2

Related Questions