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 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?

share|improve this question
2  
Unfortuantely without the base projection we can't calculate the 'per square degree' as the value will vary by large amounts depending on origin of the raster and subsequently the projection. –  om_henners Jun 25 '12 at 6:53
    
okay, i'll try to get that info. but I suppose when using GDAL or any other library, the source projection is a parameter I can chose and doesn't really affect the algorithm / steps to use. or am I wrong? –  andreash Jun 25 '12 at 12:34
    
@om A "global" raster will have its origin either at (lat,lon) = (-90,-180) or (-90,0). Andreash: the source projection is not something you choose; it is an inherent property of the data. You need to know it (or at least make an educated guess). If you guess wrong, your values could be way off, as om_henners suggests. –  whuber Jun 25 '12 at 13:58
    
Despite having posted an answer, I am now wondering if perhaps your raster represents a quantity you would prefer to express per steradian or SI squared degree (relative to the earth's center) rather than per "square degree" of latitude and longitude, as I originally interpreted it. Please let us know if this is the case, because the answer would be (quite) different! –  whuber Jun 25 '12 at 14:34
    
@whuber I know, and actually wanted to write "try to find out from the creators of the dataset which projection they used" –  andreash Jun 25 '12 at 15:31

1 Answer 1

up vote 7 down vote accepted

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. Therefore any quantity expressed as a density per square meter must be multiplied by 12.2364 E+10 * (180/pi * (sin(f+0.1) - sin(f)) / 0.1) to convert it to a density per square degree. That huge value of 12.2364 E+10 is the square of 111,194.9. (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.

share|improve this answer
3  
It seems prudent to remark that a density per square degree is a strange quantity to work with, almost without scientific meaning. For instance, something that is uniformly spread across the earth should be thought of as having constant density, whereas its density per square degree will decrease away from the equator and become practically zero at the poles. –  whuber Jun 25 '12 at 14:10

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.