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.mattiu...@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@xmax - ext@xmin)/resolution[2] > rows <- (ext@ymax - ext@ymin)/resolution[1] > > > if(inner) > { > cols <- floor(cols) > rows <- floor(rows) > } else > { > cols <- ceiling(cols) > rows <- ceiling(rows) > } > > > olim <- resolution * c(rows,cols) > ext@xmax <- ext@xmin + olim[2] > ext@ymin <- ext@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]] > > > > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
