# Easiest way to create an area raster

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 ?

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
r_sq = 6371000**2
n_lats = int(360./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