[R-sig-Geo] raster resolution

Ross Ahmed rossahmed at googlemail.com
Sun Feb 24 15:40:10 CET 2013


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.mattiuzzi at 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 <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]]
> >
> >
> >
> >
> 
> 



More information about the R-sig-Geo mailing list