1

I have a raster layer of which I would need to find the area (m2) per cell. Any suggestions for a simple approach to achieve this?

  • 1
    Project your raster layer into an appropriate coordinate reference system (CRS) that uses meters. E.g. find the correct UTM zone and use that CRS. – Jon Apr 1 at 17:17
  • 1
    @Jon that's going to change the grid basis. Probably not a good idea. If the raster is already in a cartesian system then multiply the height by the width. If the grid is lat-long then there's some calculations to be done... – Spacedman Apr 1 at 17:34
  • 1
    I was assuming the data were in lat/long. Otherwise I thought it was obvious to just look at the pixel resolution, but maybe not. – Jon Apr 1 at 17:53
  • Thanks. To clarify - with height and width, are you referring to the cell size? – A.N Apr 1 at 17:53
  • 1
    Then the unit is degrees and area of a pixel in hectares depends on where on the Earth the pixel is located. The answer by @Jon is relevant. If your raster is rather small you can calculate the length of one degree latitude and longitude around your location with some online tool like csgnetwork.com/degreelenllavcalc.html and use the result as an estimate. – user30184 Apr 1 at 19:44
0

Your question is not well constrained. If you are wondering how to find the area of a cell from a raster that is unprojected and aligned to WGS84 grid (e.g. epsg:4326), then the following Python code should help.

import numpy as np
import rasterio

fp = 'your_tiff_file.tif'
with rasterio.open(fp) as testif:
    rast = testif.read(1)
    gt = testif.transform
    pix_width = gt[0]
    ulX = gt[2]
    ulY = gt[5]
    rows = testif.height
    cols = testif.width
    lrX = ulX + gt[0] * cols
    lrY = ulY + gt[4] * rows

lats = np.linspace(ulY,lrY,rows+1)

a = 6378137
b = 6356752.3142

# Degrees to radians
lats = lats * np.pi/180

# Intermediate vars
e = np.sqrt(1-(b/a)**2)
sinlats = np.sin(lats)
zm = 1 - e * sinlats
zp = 1 + e * sinlats

# Distance between meridians
#        q = np.diff(longs)/360
q = pix_width/360

# Compute areas for each latitude in square km
areas_to_equator = np.pi * b**2 * ((2*np.arctanh(e*sinlats) / (2*e) + sinlats / (zp*zm))) / 10**6
areas_between_lats = np.diff(areas_to_equator)
areas_cells = np.abs(areas_between_lats) * q

areagrid = np.transpose(np.matlib.repmat(areas_cells,cols,1))

This snippet returns the array areagrid that contains the area in square kilometers of each pixel in your 4326 raster. There are a number of considerations "under the hood" that I am not detailing here, but I will note that I have compared this method with reprojection methods and the results have always been very similar.

You asked for a QGIS solution, so you would have to modify this a bit to get it to run in the Python console. Maybe overkill for what you want.

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.