Assign a CRS to a raster in R

Question

I have a raster that does not contain any crs. After I assign a crs, the raster still appears in another location. Am I missing something?

library(raster)
#> Loading required package: sp
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)
library(terra)
#> terra 1.7.46
library(exactextractr)
library(sp)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.1; sf_use_s2() is TRUE

r <- raster('data_bw/no_proj.tif',
                    package = 'exactextractr')
r
#> class      : RasterLayer 
#> dimensions : 10043, 8921, 89593603  (nrow, ncol, ncell)
#> resolution : 25, 25  (x, y)
#> extent     : 3388038, 3611062, 5266712, 5517788  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : no_proj.tif 
#> names      : no_proj

crs(r) <- st_crs(32632)$wkt # a WKT string

terra::writeRaster(r, "data_bw/proj.tif", filetype = "GTiff", overwrite = TRUE)

EDIT: this is the information contained in the .prj file.

PROJCS["Bessel_1841_Transverse_Mercator",GEOGCS["GCS_Bessel_1841",DATUM["D_Bessel_1841",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",3500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

Answer 1

The raster tells you its extent (the coordinate values at the corners) are:

extent     : 3388038, 3611062, 5266712, 5517788  (xmin, xmax, ymin, ymax)

Then you do this:

crs(r) <- st_crs(4326)$wkt # a WKT string

which is saying "Those coordinates are EPSG code 4326". That CRS is lat-long degrees, which is clearly wrong for those numbers.

You need to know what CRS those extent numbers are in. Likelihood is it is something like a UTM zone, in which case the numbers are metres from a 0.0 point. But without knowing this, you can't transform the CRS or otherwise georeference it.

If your TIFF is a properly-constructed GeoTIFF file then the coordinate system information should be embedded in the metadata of the .tif. If you have an associated .prj file with a coordinate system reference in it then you can read it in an assign it to the raster:

> prj = readLines("test.prj")

Raster with NA CRS:

> r
class      : RasterLayer 
dimensions : 3, 4, 12  (nrow, ncol, ncell)
resolution : 0.25, 0.3333333  (x, y)
extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 1, 12  (min, max)

Set it:

> crs(r) = prj

Now has a CRS - transverse mercator, Bessel ellipsoid etc:

> r
class      : RasterLayer 
dimensions : 3, 4, 12  (nrow, ncol, ncell)
resolution : 0.25, 0.3333333  (x, y)
extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs 
source     : memory
names      : layer 
values     : 1, 12  (min, max)