[R-sig-Geo] getting rasters to match

[Robert J. Hijmans]

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
>