# How to calculate area of 1 x 1 degree cells in a raster

I have a simple array (not geo-referenced in any way) in Matlab (which is not a language I know a huge amount about...so I can export it to other languages if needed) which has the dimensions 180 x 360 and contains data for each 1 x 1 degree area of the globe. I need to calculate the area of certain cells in this raster.

Obviously the area will vary based on longitude, but I have no idea how to calculate this area, particularly as my array is not georeferenced. I saw a reference online to the `area` function in the R `raster` package, but I'm not sure how to get my matlab data into a `Raster` object within R.

Does anyone have any suggestions as to how best to calculate these areas?

This is how you can do that with R/raster

``````library(raster)
r <- raster()  # by default 1 by 1 degree, just what you want
a <- area(r)
``````

Each cell of RasterLayer 'a' has a value representing its approximate area

To illustrate the results for one column (it is the same for all columns), as area varies by latitude, not by longitude

``````lat <- yFromRow(r, 1:nrow(r))
area <- a[,1]
plot(lat, area)
``````
• +1. The units appear to be square kilometers. BTW, `area` does not appear to be a matrix. On my system (R.2.11.1) I have to use `area <- a[seq(from=1, by=360, length.out=180)]` to extract one value per row. On the average, `area` overestimates the cell area by 0.2%. – whuber Jul 20 '12 at 12:58
• `a` is a RasterLayer. You can do `m <- as.matrix(a)` to get a matrix. I'll try your formula instead. You are using a very old version of R. It would be more useful to report the version of raster. – Robert Hijmans Jul 22 '12 at 0:26
• The 0.2% difference you found is probably because of using a different radius. Using a radius of 6378137 m, the difference is less than 0.0001% – Robert Hijmans Jul 22 '12 at 4:36
• Yes RobertH, you are right that the reason is likely the use of that radius. However, that's the wrong radius! It results in a total surface area of 511207893 km^2, which is 0.223% greater than it actually is. That's why the authalic radius of 6371007.2 meters should be used for such calculations. – whuber Jul 23 '12 at 12:06

It is a consequence of a theorem of Archimedes (c. 287-212 BCE) that for a spherical model of the earth, the area of a cell spanning longitudes l0 to l1 (l1 > l0) and latitudes f0 to f1 (f1 > f0) equals

``````(sin(f1) - sin(f0)) * (l1 - l0) * R^2
``````

where

• l0 and l1 are expressed in radians (not degrees or whatever).

• l1 - l0 is calculated modulo 2*pi (e.g., -179 - 181 = 2 degrees, not -362 degrees).

• R is the authalic Earth radius, almost exactly 6371 km.

(As a quick check, the surface area of the entire globe can be computed by letting l1 - l0 = 360 degrees = 2 Pi radians, f1 = 90 degrees, and f0 = -90 degrees. The formula gives (1 - -1) * 2 * Pi * R^2 = 4 * Pi * R^2, as is well known.)

In this raster, l1 - l0 is constantly 1 degree (0.01745329 radians) but sin(f1) - sin(f0) changes from row to row. The cell areas can therefore be computed in terms of the row indexes alone.