Convert area from square meters to degrees globally

Question

I have a global raster file with 3600x1800 pixels. I don't have any information about the projection used.

The raster contains one band with a quantity with the unit "per square meter". How can I convert this to be "per square degree", prefereably using GDAL or QGis?

Answer 1

Almost all global rasters twice as long as high, and especially those with numbers of rows and columns being nice multiples of 180, use ("unprojected") geographic coordinates. Therefore this one almost surely uses square cells of 0.1 degree. Please confirm this in the grid's metadata if possible.

The conversion from square degree to square meter depends on latitude (but not on longitude). The most accurate way to perform the conversion is to create a sequence of 0.1 by 0.1 degree polygons extending along a meridian from the equator to a pole, project those with an equal-area projection, and compute the areas in square meters. Join those areas to the grid's attribute table (using latitude as the key), then multiply the grid values by those areas and multiply the answer by 100.

We can anticipate the answers to about four significant figures by using a spherical (instead of ellipsoidal) earth model. A simple equal-area projection is cylindrical (geometrically projecting outward from the earth's axis onto an enclosing cylinder): this sends a point at (lon, lat) to (lon, sin(lat)). Thus, a square of width dl (in longitude) bounded below at a latitude f and above at a latitude f + df is projected onto a rectangle of width dl and height sin(f+df) - sin(f). Taking dl = df = 0.1 degree gives the area as

0.1 * (sin(f+0.1) - sin(f))

square degrees. (Calculus tells us that 0.1 * 0.1 * cos(f+0.05) * pi/180 is a good approximation to this, accurate almost to seven significant figures.)

To convert these degrees to meters, note that the radius of a sphere whose area is equal to that of the earth's ellipsoidal surface (an authalic sphere) is 6371 km, giving 111,194.9 meters per degree. That means a (nominal) square degree has an area of 111194.9^2 = 1.22364E+10 square meters. Therefore any quantity expressed as a density per square meter must be multiplied by 1.22364E+10 * (180/pi * (sin(f+0.1) - sin(f)) / 0.1) to convert it to a density per square degree. (Incidentally, as a double-check, this difference of sines should start around pi/180 = 0.01745 near the equator and smoothly decrease to zero at the poles; the entire expression (180/pi * (sin(f+0.1) - sin(f)) / 0.1) starts at 0.999999 at the equator where f=0 and decreases to .000873 at the pole where f=89.9. The strange quantity 180/pi is required to convert from radians to degrees.)

This calculation is straightforward to do using standard grid "map algebra" expressions (in GRASS/QGIS, ArcGIS, Manifold, Idrisi, R, or whatever): f needs to be a grid of either latitude values or row indexes and the rest is all done with pointwise ("local") operations.