raster resolution

Dear Ross,
in your mail you wrote 100 qkm and 1000 qkm. The map unit of sinusoidal is [m] as you see in the CRS.

sqrt(100)*1000 = 10000 [m] side length
and 
sqrt(1000)*1000
is the equivalent side-length of such a squared pixel in meter. The *1000 is only to convert km in m, and the sqrt() is the area to side-length conversion.
So I did not generate a raster with 1000 qm but 1000 qkm as you asked. if you want 1000 qm, resolution=sqrt(1000) 

Matteo

Ross Ahmed <rossahmed at googlemail.com> 02/24/13 9:04 AM >>>

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