Michael Dales
11/14/2022, 10:36 AMKartik Agaram
11/14/2022, 11:53 PMelevation_layer = Layer.layer_from_file('elecation.tiff')
area_layer = UniformAreaLayer('area.tiff')
validity_layer = Layer.layer_from_file('validity.tiff')
# Work out the common subsection of all these and apply it to the layers
intersection = Layer.find_intersection(elecation_layer, area_layer, validity_layer)
elevation_layer.set_window_for_intersection(intersection)
area_layer.set_window_for_intersection(intersection)
validity_layer.set_window_for_intersection(intersection)
# Work out the area where the data is valid and over 3000ft
def is_munro(data):
return numpy.where(data > 3000.0, 0.0, 1.0)
result = validity_layer * area_layer * elevation_layer.apply(is_munro)
result_band = result_gdal_dataset.GetRasterBand(1)
result.save(result_band)
I don't follow why the area_layer
is in the computation of result
. Or is it perhaps not needed in the definition of intersection
, just result
? Are the multiplications over matrices (so the ranks need to line up)? Just trying to make sure I understand the example. I don't know anything about GDAL.Konrad Hinsen
11/15/2022, 8:07 AMMichael Dales
11/15/2022, 9:46 AMresult.sum()
The number you’d get would be the square meterage of planet covered by munros. (Assuming validity layer is a layer that just covers Scotland 🙂
For the actual work we used this tiff in another set of calculations. In that work we’re reasoning about area of habitat of species, so the result in that case is a GeoTIFF per species that has a 0 value where the species isn’t, and an area of the pixel if the species is there.Konrad Hinsen
11/15/2022, 3:18 PMMichael Dales
11/16/2022, 6:32 PM