[R-sig-Geo] writeRaster and coord ref problem
Johannes Radinger
JRadinger at gmx.at
Wed Mar 7 11:29:14 CET 2012
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 at gmail.com>
> An: Johannes Radinger <JRadinger at gmx.at>
> CC: r-sig-geo at 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 at 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 at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
--
Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
More information about the R-sig-Geo
mailing list