Dear Ross,
You have 2 problems, one is that by changing the resolution you change
in most of the cases also the extent, and second you cannot generate a
regular raster (same area for all pixels) based on LAT/LON.
You will have to use to an equal area projection.
http://en.wikipedia.org/wiki/Map_projection#Equal-area
Here a function that helps you to adapt the extent:
ext=rasterExtent
resolution= rasterExtent units resolution as c(x,y) or single value for
both
inner = as the change of the resolution in most of the cases changes the
given extent you can choose to make the extent little smaller
(inner=TRUE) or little bigger.
obeyRes <- function(ext,resolution,inner=TRUE)
{
if(length(resolution)==1)
{
resolution <- rep(resolution,2)
}
cols <- (ext at xmax - ext at xmin)/resolution[2]
rows <- (ext at ymax - ext at ymin)/resolution[1]
if(inner)
{
cols <- floor(cols)
rows <- floor(rows)
} else
{
cols <- ceiling(cols)
rows <- ceiling(rows)
}
olim <- resolution * c(rows,cols)
ext at xmax <- ext at xmin + olim[2]
ext at ymin <- ext at ymax - olim[1]
return(list(extent=ext,cols=cols,rows=rows))
}
myExtent <- extent(-176.5813, 174.1103,-31.16667, 83.1)
r <- raster(ext=myExtent,crs="+proj=longlat +ellps=WGS84
+datum=WGS84")
# Sinusoidal (as example!!)
EAE <- projectExtent(r,crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=6371007.181 +units=m +no_defs")
ext100 <- obeyRes(extent(EAE),resolution=10000)
r100 <- raster(ext=ext100$extent,ncol=ext100$cols,nrows=ext100$rows,
crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181
+units=m +no_defs")
r100[] <- round(runif(ncell(r100),0,1))
# writeRaster(r100,"r100.tif")
ext1000 <- obeyRes(extent(EAE),resolution=sqrt(1000)*1000)
r1000 <-
raster(ext=ext1000$extent,ncol=ext1000$cols,nrows=ext1000$rows,
crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181
+units=m +no_defs")
r1000[] <- round(runif(ncell(r1000),0,1))
# writeRaster(r1000,"r1000.tif")
#################
I'm not 100 % sure about my calculations/considerations so take it with
care!
Matteo
Ross Ahmed 02/20/13 11:03 AM >>>
I have the following raster:
myRaster <- c(-176.5813, 174.1103, -31.16667, 83.1)
myRasterExtent <- extent(matrix(myRaster, ncol=2, nrow=2, byrow=T))
myRasterExtent <- raster(myRasterExtent, crs="+proj=longlat
+ellps=WGS84
+datum=WGS84?)
How can I change resolution to 100km^2 and 1000km^2
many thanks
Ross
[[alternative HTML version deleted]]