getting rasters to match
Advait, When I do what you did with some example data, a clear error warning message results, saying that your crs is not valid.
library(raster)
Loading required package: sp raster version 1.7-12 (9-December-2010)
futurehdr <- raster(ncol=10, nrow=10) futurehdr[] <- 1 fhdrcrs <- "+proj=lcc +datum=NAD83 +ellps=GRS80 +units=m, xmn=-2736000, xmx=2334000, ymn=-1608000, ymx=1952000" projection(futurehdr) <- fhdrcrs x <- writeRaster(futurehdr,"mask.regcm.2040-49.climatol.annsum.HDR.future.tif", overwrite=TRUE)
Geospatial Data Abstraction Library extensions to R successfully loaded Loaded GDAL runtime: GDAL 1.7.2, released 2010/04/23 Path to GDAL shared files: C:/soft/R/R-2.12.0/library/rgdal/gdal Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009 Path to PROJ.4 shared files: C:/soft/R/R-2.12.0/library/rgdal/proj Warning message: In .newCRS(value) : +proj=lcc +datum=NAD83 +ellps=GRS80 +units=m, xmn=-2736000, xmx=2334000, ymn=-1608000, ymx=1952000 is not a valid proj4 CRS string
This can be easily fixed:
crs <- "+proj=lcc +datum=NAD83 +ellps=GRS80 +units=m" projection(futurehdr) <- crs x <- writeRaster(futurehdr,"mask.regcm.2040-49.climatol.annsum.HDR.future.tif", overwrite=TRUE) x
class : RasterLayer filename : C:/Documents and Settings/rhijmans/My Documents/mask.regcm.2040-49.climatol.annsum.HDR.future.tif dimensions : 10, 10 (nrow, ncol) ncell : 100 min value : 1 max value : 1 projection : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) resolution : 36, 18 (x, y) Do note that projection(x) <- _sets_ (defines) the crs (projection) of an object. It does _not_ change it from one projection to another. You need to use projectRaster for that. so perhaps you want something like: r <- raster(vulcan.10k.dp.build.387.RES.ann1.vulcproj.img) # if necessary # projection(r) <- crs res <- projectRaster(futurehdr, r) Robert
On Sat, Dec 11, 2010 at 9:44 PM, Advait Godbole <advaitgodbole at gmail.com> wrote:
Dear List, I have a raster layer: http://dl.dropbox.com/u/4030944/vulcan.10k.dp.build.387.RES.ann1.vulcproj.img which is 356 x 507, and I want the projection and extents of this netcdf file: http://dl.dropbox.com/u/4030944/mask.regcm.2040-49.climatol.annsum.HDR.future.ncto match that. I tried doing the following: futurehdr <- raster("mask.regcm.2040-49.climatol.annsum.HDR.future.nc") fhdrcrs <- "+proj=lcc +datum=NAD83 +ellps=GRS80 +units=m, xmn=-2736000, xmx=2334000, ymn=-1608000, ymx=1952000" projection(futurehdr) <- fhdrcrs writeRaster(futurehdr,"mask.regcm.2040-49.climatol.annsum.HDR.future.tif") but it does not seem to have worked.
vulcemitres <- raster("vulcan.10k.dp.build.387.RES.ann1.vulcproj.img")
rcmtest <- raster("mask.regcm.2040-49.climatol.annsum.HDR.future.tif")
vulcemitres
class ? ? ? : RasterLayer filename ? ?: C:/Work/13Mar.research.bkup/MS_Thesis/work.data.analysis/time.series.analysis/vulcan.10k.dp.build.387.RES.ann1.vulcproj.img nrow ? ? ? ?: 356 ncol ? ? ? ?: 507 ncell ? ? ? : 180492 min value ? : 0 max value ? : 1207978 projection ?: +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 extent ? ? ?: -2735999, 2334001, -1608001, 1951999 ?(xmin, xmax, ymin, ymax) resolution ?: 10000, 10000 ?(x, y)
rcmtest
class ? ? ? : RasterLayer filename ? ?: C:/Work/13Mar.research.bkup/MS_Thesis/work.data.analysis/time.series.analysis/mask.regcm.2040-49.climatol.annsum.HDR.future.tif nrow ? ? ? ?: 356 ncol ? ? ? ?: 507 ncell ? ? ? : 180492 min value ? : ? max value ? : ? projection ?: NA extent ? ? ?: -137.3312, -62.11259, 22.12693, 57.35954 ?(xmin, xmax, ymin, ymax) resolution ?: 0.1483602, 0.09896801 ?(x, y) How can this be accomplished? Thanks, -- advait godbole ? ? ? ?[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo