[R-sig-Geo] raster resolution

[Ross Ahmed]

Dear Matteo

Your answer has been very useful - thanks. I have a query: to make a
raster in which each cell is 1000^2m, you have provided this code:

ext1000 <- obeyRes(extent(EAE),resolution=sqrt(1000)*1000)


Why have you used sqrt(1000)*1000? Why not 1000*1000?

Ross


On 20/02/2013 16:08, "Matteo Mattiuzzi" <matteo.mattiuzzi at boku.ac.at>
wrote:

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