Johannes,

Yes, that smells like something is wrong with your GDAL/PROJ4. Can you try
this:

library(rgdal)

crs <- "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
+datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0"

CRS(crs)

For me, it returns:
CRS arguments:
 +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam
+units=m +no_defs
+ellps=bessel +towgs84=606.0,23.0,413.0

Robert

On Wed, Mar 7, 2012 at 2:29 AM, Johannes Radinger <JRadinger@gmx.at> wrote:

> Hi again,
>
> I also tried your example and I can reproduce the problem in contrast to
> you. After opening a fresh session, here is what I get:
>
> > library(raster)
> Lade n"otiges Paket: sp
> raster version 1.9-70 (27-February-2012)
> > crs <- "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
> +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0"
> > r <- raster(nc=10, nr=10, crs=crs)
> > r[] <- 1:ncell(r)
> > b <- brick(r,r,r)
> > x <- writeRaster(b, 'datum.tif', overwrite=TRUE)
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.9.0, released 2011/12/29
> Path to GDAL shared files: /Users/Johannes
> Radinger/Library/R/2.14/library/rgdal/gdal
> Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
> Path to PROJ.4 shared files: /Users/Johannes
> Radinger/Library/R/2.14/library/rgdal/proj
> > x
> class       : RasterBrick
> dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
> resolution  : 36, 18  (x, y)
> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
> +ellps=bessel +towgs84=606,23,413,0,0,0,0 +units=m +no_defs
> values      : /Users/Johannes
> Radinger/Documents/Eclipse_workspace/datum.tif
> min values  : 1 1 1
> max values  : 100 100 100
>
> > sessionInfo()
> R version 2.14.2 (2012-02-29)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] rgdal_0.7-8   raster_1.9-70 sp_0.9-95     rj_1.0.0-3
>
> loaded via a namespace (and not attached):
> [1] grid_2.14.2    lattice_0.20-0 rj.gd_1.0.0-1  tools_2.14.2
>
>
>
> What might cause the problem? Is there any bug I should report? Maybe the
> problem is platform specific (I use a Mac) or has something to do with my
> version of GDAL/PROJ4?
> Any suggestions?
>
> In the meanwhile I'll try your approach with resetting the CRS:
> projection(x) <- crs
>
> /Johannes
>
>
>
> -------- Original-Nachricht --------
> > Datum: Tue, 6 Mar 2012 09:20:09 -0800
> > Von: "Robert J. Hijmans" <r.hijmans@gmail.com>
> > An: Johannes Radinger <JRadinger@gmx.at>
> > CC: r-sig-geo@r-project.org
> > Betreff: Re: [R-sig-Geo] writeRaster and coord ref problem
>
> > Johannes, it works for me:
> >
> > library(raster)
> > crs <- "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
> > +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0"
> > r <- raster(nc=10, nr=10, crs=crs)
> > r[] <- 1:ncell(r)
> > b <- brick(r,r,r)
> > x <- writeRaster(b, 'datum.tif', overwrite=TRUE)
> >
> > > x
> > class       : RasterBrick
> > dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
> > resolution  : 36, 18  (x, y)
> > extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> > coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
> > +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0
> > values      : datum.tif
> > min values  : 1 1 1
> > max values  : 100 100 100
> >
> >
> > If for some reason part of the proj.4 description is lost, you can set it
> > again with something like
> >
> > projection(x) <- crs
> >
> > > sessionInfo()
> > R version 2.14.1 (2011-12-22)
> > Platform: x86_64-pc-mingw32/x64 (64-bit)
> >
> > locale:
> > [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United
> > States.1252    LC_MONETARY=English_United States.1252
> > [4] LC_NUMERIC=C                           LC_TIME=English_United
> > States.1252
> >
> > attached base packages:
> > [1] stats     graphics  grDevices utils     datasets  methods   base
> >
> > other attached packages:
> > [1] rgdal_0.7-5   raster_1.9-71 sp_0.9-93
> >
> > loaded via a namespace (and not attached):
> > [1] grid_2.14.1    lattice_0.20-0
> >
> >
> > Best, Robert
> >
> > On Tue, Mar 6, 2012 at 6:36 AM, Johannes Radinger <JRadinger@gmx.at>
> > wrote:
> >
> > > Hi,
> > >
> > > I need to load some data from the Worldclim Database and reproject them
> > > to some other, already existing rasters. Because I don't want to use
> the
> > > reprojected rasters from the memory I want to save the brick as a
> > GeoTiff
> > > (like the other rasters) using writeRaster. But somehow writeRaster()
> > > changes the coord. ref. line and drops "+datum=Potsdam". Here what I am
> > > doing:
> > >
> > > > WC.bio <- projectRaster(getData("worldclim", var="bio", res=0.5,
> > > lon=9.3, lat=54.5),second.map)
> > > > WC.bio
> > > class       : RasterBrick
> > > dimensions  : 912, 840, 766080, 19  (nrow, ncol, ncell, nlayers)
> > > resolution  : 50, 50  (x, y)
> > > extent      : 3505300, 3547300, 6025800, 6071400  (xmin, xmax, ymin,
> > ymax)
> > > coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
> > > +datum=potsdam +units=m +no_defs +ellps=bessel
> +towgs84=606.0,23.0,413.0
> > > values      : in memory
> > > min values  :   78.0   58.0   26.0 5675.8  194.0  -24.0  214.0   52.1
> > > 28.0  151.0 ...
> > > max values  :   85   68   29 5964  210  -16  229  156   35  161 ...
> > >
> > > > writeRaster(WC.bio, "/path/to/file.tif", format="GTiff", overwrite =
> > > TRUE)
> > > class       : RasterBrick
> > > dimensions  : 912, 840, 766080, 19  (nrow, ncol, ncell, nlayers)
> > > resolution  : 50, 50  (x, y)
> > > extent      : 3505300, 3547300, 6025800, 6071400  (xmin, xmax, ymin,
> > ymax)
> > > coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0
> > > +ellps=bessel +towgs84=606,23,413,0,0,0,0 +units=m +no_defs
> > > values      : /Users/Johannes
> > > Radinger/Documents/IGB/IMPACT_Data/Treene/SDM/Predictors/WC.bio.tif
> > > min values  :   78.0   58.0   26.0 5675.8  194.0  -24.0  214.0   52.1
> > > 28.0  151.0 ...
> > > max values  :   85   68   29 5964  210  -16  229  156   35  161 ...
> > >
> > >
> > > This is a problem to my as I want to combine the other maps with the
> > > reprojected Worldclim maps into a new stack from the files. But these
> > files
> > > have different projections....
> > >
> > > Best regards,
> > > Johannes
> > >
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > R-sig-Geo@r-project.org
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
>
> --
> NEU: FreePhone 3-fach-Flat mit kostenlosem Smartphone!
> Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
>

	[[alternative HTML version deleted]]

