0

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]]

2
  • The "package" argument to raster isn't doing anything. For a second I thought this was getting a sample data file from the exactextractr package but it isn't. We can't run this code without the data file.
    – Spacedman
    Nov 27, 2023 at 11:44
  • 1
    Why are you reading the data with raster and writing with terra? The terra package is the replacement for raster so, proper workflow would be terra::rast() which results in a SpatRaster class object. Your code, as is, would not work. There is not a tif file in the exactextractr package and you are reading the data as a RasterLayer object but not coercing to a SpatRaster before using terra::writeRaster to write the new tif, I would highly recommend just dropping raster entirely and moving everything to terra. In the long run, your approach will cause numerous issues and confusion. Nov 27, 2023 at 15:57

1 Answer 1

3

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)
5
  • Thank you for your reply. I also have in the same folder a .prj and .a hdr file. Is this helpful?
    – Daniel AG
    Nov 27, 2023 at 11:48
  • this is what I see when I open 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]]
    – Daniel AG
    Nov 27, 2023 at 11:49
  • That changes everything. See edit.
    – Spacedman
    Nov 27, 2023 at 12:02
  • that's awesome, thank you! how would you, as a last step, convert the raster to epsg:32632 ?
    – Daniel AG
    Nov 27, 2023 at 12:56
  • That's a different question (to reproject the raster). Consider opening a new issue, but firstly, you should spend some time searching if that question has already being solved. Although I believe a simple google search on how to reproject raster data in R should return plenty of tutorials.
    – Nikos
    Nov 27, 2023 at 14:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.