1

The objective is to create a global geotiff raster at 30s resolution (appr 1km) with the pixel value the area of that pixel. I have multiple tools but struggle to find the easiest way. I am aware of this question and probably prefer an answer using GDAL.

QGIS:

  1. create grid,
  2. use raster calculator to obtain area per polygon
  3. rasterize featurr\es

Matlab:

  1. Areamat function

Earth Engine:

  1. ee.Image.pixelArea()

ArcMap: ?

GDAL/ Python ?

Other ?

1

Don't know if this is the easier, but using Python (Numpy) and GDAL, I think the following code works. I have calculated the grid cell size from the usual formula, and you need to provide the grid resolution in decimal degrees.

import numpy as np
import gdal
"""
Creates a global grid of a particular resolution, where the value of
each grid cell is the area of that grid cell in sqm.
"""

def do_grid (resolution=0.5):
    """Calculate the area of each grid cell for a user-provided
    grid cell resolution. Area is in square meters, but resolution
    is given in decimal degrees."""
    # Calculations needs to be in radians
    lats = np.deg2rad(np.arange(-90,90, resolution))
    r_sq = 6371000**2
    n_lats = int(360./resolution) 
    area = r_sq*np.ones(n_lats)[:, None]*np.deg2rad(resolution)*(
                np.sin(lats[1:]) - np.sin(lats[:-1]))
    return area.T

if __name__ == "__main__":
    resolution=0.5
    # Get the GeoTIFF driver
    drv = gdal.GetDriverByName("GTiff")
    # Compressed GeoTIFF file
    dst_ds = drv.Create("grid.tif", int(360./resolution),
                int(180./resolution), 
                1, gdal.GDT_Float32,
                options = [ 'COMPRESS=DEFLATE'] )
    # Projection using EPSG:4326
    wgs84='GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'
    dst_ds.SetProjection(wgs84)
    geotransform = (-180.,resolution,0,90.,0,-resolution)
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.GetRasterBand(1).WriteArray(do_grid(resolution=resolution))
    dst_ds = None
  • thank you for this answer. If I am not mistaken, this will only work with rasters smaller than your memory. – RutgerH Mar 23 '17 at 7:31
  • Well, you can just create chunks of the raster and write them out. Or you could create a coarser resolution raster and then use gdalwarp it to the desired resolution. This second option is easier I think. – Jose Mar 24 '17 at 9:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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