3

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?

A.N
  • 181
  • 6
  • 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 01 '19 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 01 '19 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 01 '19 at 17:53
  • Thanks. To clarify - with height and width, are you referring to the cell size? – A.N Apr 01 '19 at 17:53
  • The cell size means the same as pixel size. You can find that with QGIS from Layer-Properties-Information. – user30184 Apr 01 '19 at 18:33
  • 1
    Okay, thanks. My pixel size is 0.000833333, -0.000833333. As I understand, this is not in the unit of m2? I have a raster containing values of energy demand per hectare per cell. I would like to convert this to energy demand per cell and thus need to identify the area of the cell. – A.N Apr 01 '19 at 18:41
  • 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 http://www.csgnetwork.com/degreelenllavcalc.html and use the result as an estimate. – user30184 Apr 01 '19 at 19:44
  • Okat I see, many thank for this it should be very helpful. Although I have one last question; In the pixel size displayed under layer information, which value corresponds to longitude vs latitude? – A.N Apr 02 '19 at 08:09

2 Answers2

5

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.

Jon
  • 2,894
  • 2
  • 19
  • 35
0

It is possible using the SAGA raster calculator per Getting latitude / longitude in the QGIS Raster Calculator and Creating direction raster in QGIS?

Since you know your cell size is 0.000833333, -0.000833333 degrees (per comment), (or 3 arc-seconds/cell) then, you can use the xpos() and ypos() to convert the degrees into an area.

Assuming a minute of latitude is a standard nautical mile of 1852 meters on a spherical earth:

  area_m2 = 0.000833333*60*1852 * 0.000833333*60*1852* cos(ypos())`

or

 area_m2 = (0.000833333*60*1852)^2*cos(ypos()*pi()/180)

In the Saga raster calculator, (under Processing/Saga/Raster Calculus/Raster Calculator), choose your layer and use a formula like:

(0.000833333*60*1852)^2*cos(ypos()*pi()/180)

Screenshot of QGIS SAGA Raster Calculator config screen

This should produce a new raster layer conformant with your input layer, with each cell holding the square meters of a 3-arc-second cell at that latitude (about 6830m2 in my area, decreasing northward)

Dave X
  • 1,639
  • 13
  • 25