Calculate cell/pixel area with Rasterio in Python

Question

The Situation

I have a tile spanning 10° north-south (n-s) and 10° east-west (e-w), e.g. from 10-20W and 10-20N. This tile contains multiple pixels, e.g. with a resolution of n-s and e-w of 100m at the equator (~0.0008°). This raster contains information of biomass in tons/hectare (for every pixel).

These 10x10° tiles span the whole world and are currently in EPSG:4326 - WGS 84.

In the end, I would like the absolute biomass in every pixel, and therefore need to multiply the tons/hectare information with the corresponding pixel size.

I am currently using this Determine area of cell in raster (QGIS) to create a new raster, and then do raster multiplication with the tons/hectare raster.

The Question

How can I calculate the area of every 100x100m pixel in Rasterio? I am looking for an approach that works world-wide. If possible also without creating a new raster, but rather multiplying the tons/hectare raster "inplace" with its pixel area.

I would already be happy if the current code I am using could be "translated" to Rasterio, as I am a bit confused between GDAL's raster.GetGeoTransform() and Rasterio's raster.transform.

Answer 1

The correspondence between GDAL's raster.GetGeoTransform() and Rasterio's raster.transform is easy (look at Python affine transforms and how to take cell size from raster using python or GDAL or RASTERIO)

Pixel size with rasterio

import rasterio as rio
ww = rio.open("wgs84_raster.tif"
gt = ww.transform
gt
Affine(0.0004977676400361997, 0.0, 5.096854802788522,
   0.0, -0.00031309240744181374, 50.123699567610515)
# this is an affine matrix
print(gt)
| 0.00, 0.00, 5.10|
| 0.00,-0.00, 50.12|
| 0.00, 0.00, 1.00|
pixelSizeX = gt[0]
pixelSizeY =-gt[4]
print(pixelSizeX) # in degrees
0.0004977676400361997
print(pixelSizeY) # in degrees
0.00031309240744181374

The resulting area is in square degrees. If you want square kilometers, you need to reproject the raster

Convert the rasterio transform matrix to Gdal geotransform

 print(gt.transform.to_gdal())
 (5.096854802788522, 0.0004977676400361997, 0.0, 50.123699567610515, 0.0, -0.00031309240744181374)

Pixel size with gdal

from osgeo import gdal
raster= gdal.Open("wgs84_raster.tif")
got = raster.GetGeoTransform()
print(got) 
(5.096854802788522, 0.0004977676400361997, 0.0, 50.123699567610515, 0.0, -0.00031309240744181374)

pixelSizeX = got[1]
pixelSizeY =-got[5]
print(pixelSizeX) # in degrees
0.0004977676400361997
print(pixelSizeY) # in degrees
0.00031309240744181374

The resulting area is in square degrees

Convert gdal geotransform to rasterio transform matrix

af = rio.Affine.from_gdal(*got)
af
Affine(0.0004977676400361997, 0.0, 5.096854802788522,
   0.0, -0.00031309240744181374, 50.123699567610515)

You can also use directly the affine module (included in rasterio)

from affine import Affine
af = Affine.from_gdal(*got)

But the resulting area is in square degrees and not square meters