Hi all, I am trying to calculate the area of urban land cover only for cells in which 100% of the cell is urban cover.
I am calculating the area as the product of the resolution times the number of cells, but it is a bit clumsy. My questions are:
Is my calculation correct?
Is there a better way to calculate area?
How can I calculate the area of cells in which 100% of the cell is urban cover?
Thank you for the help!
using Rasters
using RasterDataSources
using DataFrames
# Set path to download raster
ENV["RASTERDATASOURCES_PATH"] = "C:/path/"
# Get urban land cover
urban = Raster(EarthEnv{LandCover}, 9, missingval = -9999)
# Calculate resolution and area
urban_df = DataFrame(values(urban))
res = (urban_df.X[end] - urban_df.X[1])/length(urban[:, 1, 1])
prod = res * res
ncells = length(urban[:, 1, 1]) * length(urban[1, :, 1])
area = ncells * prod
There is nothing wrong, or inaccurate, to compute total area by summing the area of each pixel that fulfills a certain condition (here, being a land cover pixel). The land cover files come, if I’m not mistaken, in GeoTIFF. So the pixel resolution is part of the file’s metadata and you don’t need to compute it yourself.
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.
Right, I was assuming the GeoTIFF were in UTM (they often are) but if they are in geographical coordinates then the size of each pixel is constant along a scan line but changes from scan-line to scan-line in function of cosinus(lat).
That is, for each lat strip the pixel’s area ~=
EarthRadius * dlat * pi /180 * EarthRadius * dlon * pi /180 * cos(lat * pi/180)