Calculate area of a raster

First, looking at your code, there are a minor few errors. You are taking the first and last values as the edges of the raster, rather than the bounds - this means you are missing one pixel at the end - these are Intervals and the key values are the start of each interval. You can get the bounds of a raster by doing bounds(urban, (X, Y)).

And if a dimension has Regular spacing (as these and anything a GeoTIFF can represent has) you can get the resolution directly from the dimensions using the Base method step, and taking the absolute value with abs (the Y axis often has a negative step size in spatial data).

So:

using Rasters
using RasterDataSources
urban = Raster(EarthEnv{LandCover}, 9)[Band(1)]
res = abs.(step.(dims(urban, (X, Y))))
pixel_area = prod(res)
# Create a Bool array with anything over 50 as true (or choose another cutoff)
urban_pixels = urban .> 50
# Sum the urban pixels and multiply by the pixel area
sum(urban_pixels) * pixel_area

But, actually, you cant get an meaningful area from these numbers because EarthEnv is projected in arc seconds! It’s harder to get area in meters because your Y values will vary per pixel, and the X values will only apply at the equator. (Calculating small areas near the equator we can often ignore this, but this is the whole world, and you will greatly over state the relative size of urban areas in temperate regions.)

You will need to convert to a projection in meters, and calculate the size of every pixel using the projection. Maybe GDAL has some tools to do this.

1 Like