5

I would like to resample a raster from a high resolution to a low resolution (with different extent) in a defined grid cell. Is there a way of using an existing raster file as input for the snapping?

In the raster package, aggregate and resample seem to be adequate but I can't find how to do it.

Wraf
  • 747
  • 2
  • 10
  • 24
  • Do you need to interpolate onto a different grid, or what? All `raster` files, so far as I know, define data on a uniform rectangular grid, so "using an existing raster file" would just mean "aggregating or interpolating from MxN to LxK grid." – Carl Witthoft Oct 08 '13 at 11:25
  • The question lacked clarity and had no example. – IRTFM Oct 08 '13 at 15:39

3 Answers3

7

You can use projectRaster for this if you have a raster in one projection and resolution and you need output in a different particular resolution and projetion.

The from argument is your high resolution raster and the to argument is your low res raster. Make sure you choose the correct method for aggregation (i.e. bilinear for continuous data and ngb (nearest neighbour) for categorical data.

require( raster )

#  Projection info
proj1 <- CRS("+proj=laea +lon_0=20 +lat_0=5 +ellps=sphere +unit=km +to_meter=1e3")
proj2 <-  CRS("+proj=longlat +datum=WGS84 +ellps=WGS84")
#  High res raster
r1km <- raster( nrows = 1515 , ncols = 2300 , xmn = -4000 , xmx = -1700 , ymn = -15 , ymx = 1500 , crs = proj1 )

#  Low res raster
r5km <- raster( nrows = 303 , ncols = 460 , xmn = -20 , xmx = -5 , ymn = 4 , ymx = 15 , crs = proj2 )

#  Set some values in high res raster
pts <- rasterToPoints(r1km)
values( r1km ) <-  0.01*pts[,1] + sin(0.02*pi*pts[,2])

#  Reproject using the attributes of the low res raster for output
out <- projectRaster( from = r1km , to = r5km , method = "bilinear" )

#  Plot - extent of second raster doesn't fully cover first so some data is missing
par( mfrow = c(1,2) )
plot( r1km )
plot( out )

enter image description here

If your input and output data are the same except in resolution you can use aggregate...

#  If same extent and resolution require use aggregate
r1 <- raster(system.file("external/rlogo.grd", package="raster"))
r5 <- aggregate( r1 , fact = 5 , method = "bilinear" )
par( mfrow = c(1,2) )
plot( r1 )
plot( r5 )

enter image description here

Simon O'Hanlon
  • 58,647
  • 14
  • 142
  • 184
  • 1
    Thanks but as my two files are in the same projection, the project does not work. About the aggregate, I cannot specify a file directly in order to snap to the grid, because I don't want to use a "fact" argument but a file? – Wraf Oct 08 '13 at 09:41
  • @Wraf you should really provide a reproducible example showing what you have done then. I don't think I expend any more effort on this guessing at your requirements. Please read [**how to make a great reproducible example**](http://stackoverflow.com/q/5963269/1478381) and update your question accordingly! – Simon O'Hanlon Oct 08 '13 at 09:57
  • Thanks. Indeed, not easy to find an answer when not reproducible example is given. But the dataset used isquite big and making such example is not that easy. – Wraf Oct 08 '13 at 10:41
  • @Wraf Check out how I made an example raster in my post from scratch (or from file). Make some data up to illustrate the problem. – Simon O'Hanlon Oct 08 '13 at 10:44
1

You can launch an external command with system and call a gdal_translate or gdal_warp command. This requires of course the installation of gdal utilities

WAF
  • 1,141
  • 20
  • 44
  • Thank for the idea but using gdal_translate does not allow to control the interpolation method. I found the solution using gdalwarp – Wraf Oct 08 '13 at 10:36
1

This solution works :

system(paste("gdalwarp"
,paste(dir_path,"fileHR.tif",sep="")
,paste(dir_path,"fileLR.tif",sep=""),sep=" "))

where dir_path is the directory where you file are stored, fileHR.tif is the High resolution file, fileLR.tif is the low resoltion file.

Wraf
  • 747
  • 2
  • 10
  • 24