[R-sig-Geo] writeRaster and coord ref problem
Johannes Radinger
JRadinger at gmx.at
Fri Mar 9 11:48:25 CET 2012
Hi Robert,
thank you for helping me sourcing the problem.
Here what I get (there is not problem with CRS(crs)):
from a totally fresh session:
> library(rgdal)
Lade n"otiges Paket: sp
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
> 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)
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
> 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 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
My rgdal package for Mac OSX was installed as mentioned on the CRAN-website (http://cran.r-project.org/web/packages/rgdal/index.html)
anything else I can/should test...?
/johannes
PS: @Roger: I set you CC as you are the maintainer of the rgdal-package and might be interested in that Problem.
-------- Original-Nachricht --------
> Datum: Thu, 8 Mar 2012 22:55:34 -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,
>
> 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 at 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 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
> > > >
> >
> > --
> > NEU: FreePhone 3-fach-Flat mit kostenlosem Smartphone!
> > Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
> >
--
Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
More information about the R-sig-Geo
mailing list