[R-sig-Geo] raster resolution

[Matteo Mattiuzzi]

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]]