10

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?

Improve this question

2 Answers 2

Reset to default
17

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.

Improve this answer
9

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)
Improve this answer
4
  • +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
    Commented Jul 20, 2012 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
    Commented Jul 22, 2012 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
    Commented Jul 22, 2012 at 4:36
  • 1
    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
    Commented Jul 23, 2012 at 12:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.