0

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?

  • GeoTIFF files, which is what you get if you writeRaster(foo, "foo.tif") do store projection info, – Spacedman Aug 14 '18 at 13:40
  • @Spacedman- the written rasters in tif format does not have the datum information i.e +datum=NAD83 when read back into R or when viewed in QGIS or Arc (comes up as D_unknown) – tg110 Aug 14 '18 at 13:49
  • The information about the datum is written to the TIFF, and read back in, its just the name of the datum that isn't recognised. That doesn't matter much. – Spacedman Aug 14 '18 at 14:06
  • That is hard to believe! And very misdirecting since there is no mention of the datum in the projection string. So if I were to pass this product on to somebody else, they would think it is wrong because the datum info does not come up in the string when read into R or QGIS and comes up as unknown in Arc. Also as mentioned in the post linked in the question, including +datum and writing out the shps works, so why cannot it work with rasters? – tg110 Aug 14 '18 at 14:32
  • The datum info DOES come up in the string. Its the +ellps and the +towgs84 parameters. The NAD83 datum is defined here: epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::4269. – Spacedman Aug 14 '18 at 14:40
2

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 
| improve this answer | |

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.