Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
add comment

2 Answers 2

up vote 3 down vote accepted

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)
share|improve this answer
    
+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. –  RobertH 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% –  RobertH 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
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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