Dear Matteo Sorry for confusing metres and kilometres. If I understand correctly:
sqrt(1000)*1000 = each cell in raster is 1000^2km Therefore: sqrt(50)*50 = each cell in raster is 50^2km sqrt(300)*300 = each cell in raster is 300^2km sqrt(800)*800 = each cell in raster is 800^2km Could you let me know if these are correct? Ross On 24 Feb 2013, at 10:16, "Matteo Mattiuzzi" <matteo.mattiu...@boku.ac.at> wrote: > 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 <rossah...@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.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