1

I have a normal RasterLayer in R and I would like to write it to a file and visualize it in QGIS. When I look at in in R it says:

class      : RasterLayer 
dimensions : 3035, 1469, 4458415  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : 386042.8, 404405.3, 6293348, 6331286  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 2.220446e-16, 0.9845922  (min, max)

The I do a writeRaster(<rasterlayer>, <file.tif>, format = "GTiff", overwrite=TRUE). When I check the coordinate reference with gdalsrsinfo on the command line I get this information:

Driver: GTiff/GeoTIFF
Files: no_filter.tif
Size is 1469, 3035
Coordinate System is:
BOUNDCRS[
    SOURCECRS[
        PROJCRS["WGS 84 / UTM zone 19S",
            BASEGEOGCRS["WGS 84",
                DATUM["World Geodetic System 1984",
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4326]],
            CONVERSION["UTM zone 19S",

However, when I drag and drop it to QGIS I get the Unknown CRS output. Is there anything special that should be considererd when writing raster-files or something like a "best practice"?

UPDATE SO when I just print the raster with raster_objects[[1]] for example I get:

> raster_objects[[1]]
class      : RasterLayer 
dimensions : 3035, 1469, 4458415  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : 386042.8, 404405.3, 6293348, 6331286  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 2.220446e-16, 0.9845922  (min, max)

When I do a rasterWrite(raster_objects[[1]], "test/test.sdat", format="SAGA"), I get:

class      : RasterLayer 
dimensions : 3035, 1469, 4458415  (nrow, ncol, ncell)
resolution : 12.5, 12.5  (x, y)
extent     : 386042.8, 404405.3, 6293348, 6331286  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /home/robin/geodata/hiwijason/workPkg_LslPrd_DEMsmoothing/maipo_basin/data/saga/raster_predictions/test.sgrd 
names      : layer 
values     : 2.220446e-16, 0.9845922  (min, max)

And I have no clue why the crs-line is saying NA in the writing-operation

UPDATE Maybe this has to do with the raster beeing in SAGA-Format (.sdat)?! Because when I drag and drop the .sdat-file into qgis it gets georeferenced directly. However, when I import it into R and write it out again with:

r = raster("...sdat")
writeRaster(r, "...tif", format="GTiff")

I lose the crs somewhere on the way. Is there anything I'm doing completely wrong? Would readGdal maybe offer help?

3
  • 1
    I cannot reproduce this problem (on Windows10; raster 3.1-5, rgdal 1.5-10, tested in QGIS 3.12.3). Potentially PROJ related issue? Migration to PROJ6/GDAL3 – Kazuhito Jun 20 '20 at 17:21
  • 1
    Is that the full Coordinate System output from gdalsrsinfo? Because it looks horribly truncated - it ends with a dangling comma and lots of closing brackets are missing. – Spacedman Jun 21 '20 at 10:17
  • No no it's not the full output. I just wanted to show, that there is some output, but in Qgis there is no spatial reference. – Lenn Jun 21 '20 at 10:53
2

I ran into the same issue. The following method seems to have worked - the coordinate system showed in QGIS as well.

theCRS <- sp::CRS('+init=epsg:4283')    # GDA94 (Australia)
crs(theRaster) <- theCRS    
writeRaster(theRaster, ...)

The drawback is that you need to know the EPSG code.

If you assign the CRS (this way) to the source data before doing a calculation then the output also seems to have a recognized crs.

1

I faced this type of issue when importing and write ESA SNAP .dim images in R. Then after some research I tried to update the crs of the raster from here. It works for me.

try this for your case:

crs(raster_objects)<-CRS("+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs")

then export your image. you can find the proj4string it here.

3
  • Thank you very much:)! I tried it, but didn't change anything – Lenn Jun 27 '20 at 13:25
  • I think something is wrong with my data. I tried it with different rasters and it works just normal. Is there any option to check if the metadata of my raster is "valid" ? – Lenn Jun 27 '20 at 13:32
  • can you share this image with me? or a link? you can find me Here. – Subhadip Datta Jun 27 '20 at 14:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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