Resolution values for Rasters in R

Question

I was just wondering if anyone has a conversion table for the resolution of rasters in R. I would like to know what numbers like these refer to in meters:

resolution : 0.08333333, 0.08333333

resolution : 0.009398496, 0.009398496
resolution : 0.002349, 0.002349 = 250m (I think)

I would really like to know what resolution to set a raster object to make cell sizes 1km^2. I am using rasters that span the country of Australia.

Thanks in advance everyone.

Cheers, Adam

Answer 1

It all depends on the units of your raster, and that depends on the projection. Rasters might not even be square grids in metres - they might be square in degrees which aren't square in metres!

1 degree longitude at the equator is 1/360 of the earth's circumference. Near the north pole 1 degree is a much smaller distance, and at the pole its pretty much zero. Degrees of latitude however are constant.

You could take the corner points of your raster, convert them to lat-long coordinates if not already, and then work out the great-circle distance between them (there's an rdist function somewhere that does this I recall). However this won't work if your raster spans the whole globe, since then your NW corner and your NE corner are at the same point... Ummm. Anyway, the answer is... 42.

If you want to make 1km rasters of Australia then.... you need a coordinate system of Australia in kilometres. In the UK we have a system called the OSGB National Grid, which is close enough to a metric grid. Australia might be trickier because it is slightly bigger than the UK... So Australia seems to have a few grid systems. See here:

http://www.spatialreference.org/ref/?search=AGD84

So you might want to use the system that is in the middle of the country to avoid the worst distortions, then work out the bounds of Australia in lat-long, convert to epsg:20353 and create a raster based on that:

In lat-long I reckon Australia is roughly:

> xtll
         [,1]       [,2]
[1,] 112.5162 -43.906900
[2,] 155.8425  -7.585619

make this into a SpatialPoints object:

> xtll=SpatialPoints(xtll,CRS("+init=epsg:4326"))

convert to that AGD84 in the middle of the country:

> spTransform(xtll,CRS("+init=epsg:20353"))
SpatialPoints:
     coords.x1 coords.x2
[1,]  -1306200   4886041
[2,]   2849956   9103124

Make a raster extent object rounded to km:

> ext = extent(-1306000,2850000,4886000,9103000)

How many rows and columns do we need?

> length(-1306:2850)
[1] 4157
> length(4886:9103)
[1] 4218

Create a raster:

> r = raster(ext,ncol=4156,nrow=4217,crs="+init=epsg:20353")
> r
class       : RasterLayer 
dimensions  : 4217, 4156, 17525852  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : -1306000, 2850000, 4886000, 9103000  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:20353 +proj=utm +zone=53 +south +ellps=aust_SA +units=m +no_defs 
values      : none

Note that the ncol and nrow values are one less than the values from the bounds - this would be a fencepost error to put those values in.

See how my resolution is 1000? This is a 1km grid. The problem is that this is possibly going to be a bit distorted on the coasts. You could work out how distorted by converting to lat-long (epsg:4326), then to the proper AGD zone for points on the coast, and seeing how different they are. They might be very close, except for an offset.

Anyway, nuff said.

Answer 2

Finding the meta-data that gives meaning to your raster can be a bit of a challenge. I have spent lots of time hunting for this. If the raster was published by a government agency, then I would hope that this information is posted somewhere prominently.

The good news is that once you know the projection used on the various rasters, you can convert them to a common projection using projectRaster() in the raster package. You need to find the proj.4 string describing the original and the desired projections in each case. You can get this from: http://www.spatialreference.org.

When you know your projection, then the resolution information you seek will have meaning.