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?

whuber
  • 69,783
  • 15
  • 186
  • 281
robintw
  • 4,006
  • 11
  • 40
  • 60

2 Answers2

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.

Community
  • 1
whuber
  • 69,783
  • 15
  • 186
  • 281
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)
Robert Hijmans
  • 10,683
  • 25
  • 35
  • +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 00: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 04: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 Jul 23 '12 at 12:06