writeRaster() not writing datum information

Question

I am trying to reproject a raster to the USGS version of the Albers Equal Area Conic projection. Following this post, which talks about projecting vectors by adding the +datum=NAD83 and writing it out, I did the same for the raster and wrote it out using writeRaster().

However, the resulting raster still has no NAD83 datum when viewed in Arc or QGIS or when read back into R (i.e. the +datum=NAD83 doe snot appear in the projection string). How do I project the raster to the USGS Albers Equal Area Conic projection with the NAD83 datum and write the raster out?

Also, I am writing out the output in tif format, which I don't think can store projection information. Should I be writing out the raster in another format that can store projection information? Which one?

Answer 1

I have been able to reproduce this on my machine and have found a solution. The "data" file is quite large so I don't post it here but as an example, see below. There is some raster file data which is loaded into the workspace. I perform some calculations, and save a copy. "data" has the crs given below.

> writeRaster(data, "~/R_work/Mackenzie/2015_06_binary.tif", 
+             "GTiff", datatype = "INT1U", overwrite = T, options = c("COMPRESS=LZW", "TFW=YES"))
> data
class       : RasterLayer 
dimensions  : 24099, 13793, 332397507  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : -15267810, -14854020, 10270650, 10993620  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
names       : layer 
values      : 0, 1  (min, max)

After reading the dataset back in and looking at the crs, it appears different. This difference doesn't allow it to be read into ArcGIS or Python.

> data <- raster("~/R_work/Mackenzie/2015_06_binary.tif")
> crs(data)
CRS arguments:
 +proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs 

However, if we use the CRS function from sp (not crs from raster) this seems to be fixed. Note that the crs of data3 is missing the EPSG - this is fine, as the other properties are preserved. Observe:

data2 <- data
crs(data2) <- CRS('+init=EPSG:3857')
> crs(data2)

CRS arguments:
 +init=EPSG:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +no_defs 

writeRaster(data2, "C:/Users/lvulis/Documents/deltas/R_work/Mackenzie/2015_06_binary2.tif", 
            "GTiff", datatype = "INT1U", overwrite = T, options = c("COMPRESS=LZW", "TFW=YES"))
data3 <- raster("C:/Users/lvulis/Documents/deltas/R_work/Mackenzie/2015_06_binary2.tif")
> crs(data3)
CRS arguments:
 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m
+nadgrids=@null +no_defs