projecExtent
returns a RasterLayer with a projected extent, but without any values. This RasterLayer can then
be adjusted (e.g. by setting its resolution) and used as a template 'to'
Raster object in projectRaster
.
projectRaster
computes values for the cells of the new Raster object.projectRaster(from, to, res, crs, method="bilinear", filename="", ...)
projectExtent(object, crs)
crs
argument, and, optionally, a res
argument, but do not provide a to
argument.
2) Create a template Raster with the CRS you want to project to. You can use an existing object, or use projectExtent
for this
or an existing Raster* object. Also set the number of rows and columns (or the resolution), and perhaps adjust the extent. The resolution of the output raster should normally be similar to that of the input raster. Then use that object as from
argument to project the input Raster to.
This is the preferred method because you have most control. For example you can assure that the resulting Raster object lines up with other Raster objects.
Projection is performed using the PROJ.4 library accesed through the rgdal package.
One of the best places to find PROJ.4 coordinate reference system descriptions is projInfo('proj')
, projInfo('ellps')
, and projInfo('datum')
for valid PROJ.4 values.
The following additional arguments can be passed, to replace default values for this function
overwrite
Logical. If TRUE
, "filename" will be overwritten if it exists
format
Character. Output file type. See writeRaster
datatype
Character. Output data type. See dataType
progress
Character. "text", "window", or "" (the default, no progress bar)
}resample
CRS-class
, projInfo
, spTransform
# create a new (not projected) RasterLayer with cellnumbers as values
r <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, ncols=40, nrows=40)
r <- setValues(r, 1:ncell(r))
projection(r)
# proj.4 projection description
newproj <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
# we need to have rgdal installed for this
if (require(rgdal)) {
#simplest approach
pr1 <- projectRaster(r, crs=newproj)
# alternatively also set the resolution
pr2 <- projectRaster(r, crs=newproj, res=20000)
# inverse projection, back to the properties of 'r'
inv <- projectRaster(pr2, r)
# to have more control, provide an existing Raster object, here we create one
# using projectExtent (no values are transferred)
pr3 <- projectExtent(r, newproj)
# Adjust the cell size
res(pr3) <- 200000
# now project
pr3 <- projectRaster(r, pr3)
# using a higher resolution
res(pr1) <- 10000
pr <- projectRaster(r, pr1, method='bilinear')
inv <- projectRaster(pr, r, method='bilinear')
dif <- r - inv
# small difference
plot(dif)
# Meuse data to long/lat
meuse <- raster(system.file("external/test.grd", package="raster"))
meusell <- projectRaster(meuse, res=1/1200, crs="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
}
Run the code above in your browser using DataLab