Oliver,

I think starting your script with this should do it:

library(rgdal)
set_OVERRIDE_PROJ_DATUM_WITH_TOWGS84(FALSE)

Robert



On Tue, Nov 6, 2012 at 9:41 PM, Oliver Soong <osoong+r@gmail.com> wrote:

> So I'm also using the raster package, and the datum problems with
> rgdal are inherited.  For those who haven't been following this
> thread, there is a workaround as of rgdal 0.7.22 and setting
> OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = FALSE.  Is there any way to have
> raster respect this setting?  For small rasters, I know I can use
> readGDAL and coerce the SpatialGridDataFrame to a RasterLayer, but
> this requires being able to load the entire raster into memory.
>
> R 2.15.1 32-bit, rgdal 0.7.22, raster 2.0.12, Windows XP/7.
>
> Test code, translated from the rgdal example:
>
> ext <- extent(-2100050, -2099950, 1199950, 1200050)
> p4s <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
> +y_0=0 +datum=NAD83 +units=m +no_defs"
> img1 <- raster(nrows = 1, ncols = 1, crs = p4s, ext = ext)
> values(img1) <- 0
> img1.file <- file.path(tempdir(), "img1.tif")
> writeRaster(img1, img1.file, "GTiff", overwrite = TRUE)
> img2 <- raster(img1.file)
> img3 <- readGDAL(img1.file, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = FALSE)
> projection(img1)
> projection(img2)
> projection(img3)
>
> Thanks,
> Oliver
>
>
> On Fri, Nov 2, 2012 at 1:03 AM, Oliver Soong <osoong+r@gmail.com> wrote:
> > R 2.15.1 32-bit, rgdal 0.7.20, Windows 7.
> >
> > grid <- GridTopology(c(-2100000, 1200000), c(100, 100), c(1, 1))
> > p4s <- CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96
> > +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80
> > +towgs84=0,0,0")
> > img1 <- SpatialGridDataFrame(grid, data.frame(band1 = 1), p4s)
> > img1.file <- file.path(tempdir(), "img1.tif")
> > writeGDAL(img1, img1.file)
> > img2 <- readGDAL(img1.file)
> > img2.file <- file.path(tempdir(), "img2.tif")
> > writeGDAL(img2, img2.file)
> > img1@proj4string
> > img2@proj4string
> >
> > For me, img1@proj4string has +datum=NAD83 and img2@proj4string does
> > not.  Not surprisingly, if I look at both files in Arc, img1 has a
> > defined datum and img2 does not.
> >
> > Am I doing anything wrong?
> >
> > Oliver
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

