Skip to content

Instantly share code, notes, and snippets.

@NyakudyaA
Last active March 24, 2023 04:29
Show Gist options
  • Save NyakudyaA/b4640ec9d2b5f43fa456083b61cfd12f to your computer and use it in GitHub Desktop.
Save NyakudyaA/b4640ec9d2b5f43fa456083b61cfd12f to your computer and use it in GitHub Desktop.
PyQGIS Implementation of generating raster statistics
import numpy as np
from PyQt5.QtCore import QFileInfo
import gdal
import osr
from PyQt5.QtGui import QColor
import processing
__author__ = 'admire'
from qgis.core import QgsRasterShader, QgsColorRampShader, QgsSingleBandPseudoColorRenderer, \
QgsProject, QgsRasterLayer, \
QgsRasterBandStats, QgsCoordinateReferenceSystem
class ReportClass(object):
def __init__(self):
self.colour_ramp = {}
def mini_style(self, layer):
raster_shader = QgsRasterShader()
color_ramp = QgsColorRampShader()
color_ramp.setColorRampItemList(self.colour_ramp['ramp'])
color_ramp.setColorRampType(QgsColorRampShader.Interpolated)
raster_shader.setRasterShaderFunction(color_ramp)
my_pseudo_renderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(),
layer.type(),
raster_shader)
layer.setRenderer(my_pseudo_renderer)
layer.triggerRepaint()
QgsProject.instance().addMapLayer(layer, True)
if layer.isValid() is True:
print("Layer was loaded successfully")
else:
print("Unable to read basename and file path - Your string is probably invalid")
def get_area(self, raster_path):
raster_layer = QgsRasterLayer(raster_path)
colour_ranges = self.colour_ramp['ranges']
colours = self.colour_ramp['colours']
area_per_colour = {}
# open the raster
ds = gdal.Open(raster_path)
if ds is None:
gdal.sys.exit("ERROR: can't open raster")
# get cell sizes
# TODO: Add this piece of code correctly
pixel_size_x = raster_layer.rasterUnitsPerPixelX()
pixel_size_y = raster_layer.rasterUnitsPerPixelY()
hectare_per_cell = pixel_size_x * pixel_size_y / 10000
# get raster bands
# hectare_per_cell = 225
raster = ds.ReadAsArray()
# process the raster
for (bottom, top, colour) in zip([-1] + colour_ranges[:-1], colour_ranges, colours):
# TODO: break up
# TODO: convert to method
# TODO: test
raster = np.where(raster != -9999, raster, np.nan)
# print top
# print np.nanmin(raster)
raster_bounds = np.where(top >= raster, raster, np.nan)
# print np.nansum(raster_bounds)
lower_bounds = np.where(raster_bounds > bottom, 1, np.nan)
classified_area = np.nansum(lower_bounds)
# print classified_area
if np.isnan(classified_area):
classified_area = 0
area_per_colour[colour] = round(hectare_per_cell * classified_area, 2)
# print area_per_colour
return area_per_colour
def get_raster_area(self, raster_path):
raster_layer = QgsRasterLayer(raster_path)
ds = gdal.Open(raster_path)
if ds is None:
gdal.sys.exit("ERROR: can't open raster")
# get cell sizes
raster_band = ds.GetRasterBand(1)
no_data_value = raster_band.GetNoDataValue()
if no_data_value is None:
no_data_value = -9999
else:
no_data_value = no_data_value
# TODO: Add this piece of code correctly
pixel_size_x = raster_layer.rasterUnitsPerPixelX()
pixel_size_y = raster_layer.rasterUnitsPerPixelY()
hectare_per_cell = pixel_size_x * pixel_size_y / 10000
raster = ds.ReadAsArray()
all_pixels_1 = np.where(raster != no_data_value, 1, np.nan)
total_pixel_area = np.nansum(all_pixels_1)
total_area_raster = total_pixel_area * hectare_per_cell
return total_area_raster
# function to style raster layers dynamically
def style_raster(self, raster_path):
file_info = QFileInfo(raster_path)
path = file_info.filePath()
basename = file_info.baseName()
layer = QgsRasterLayer(path, basename)
provider = layer.dataProvider()
extent = layer.extent()
ver = provider.hasStatistics(1, QgsRasterBandStats.All)
stats = provider.bandStatistics(1, QgsRasterBandStats.All, extent, 0)
self.colour_ramp = {'colours': [], 'ranges': [], 'ramp': []}
if ver is not False:
print("minimumValue = ", stats.minimumValue)
print("maximumValue = ", stats.maximumValue)
if stats.minimumValue < 0:
minimum = 0
else:
minimum = int(stats.minimumValue)
maximum = int(stats.maximumValue)
raster_range = maximum - minimum
total_area_raster = self.get_raster_area(raster_path)
print(total_area_raster)
print(" 6 classes generated")
colours = ['#33a02c', '#80c02e', '#cee130', '#f1c62c', '#ea7024',
'#e31a1c']
add = raster_range // (len(colours) - 1)
int_1 = int(minimum + add)
int_2 = int(int_1 + add)
int_3 = int(int_2 + add)
int_4 = int(int_3 + add)
value_list = [minimum, int_1, int_2, int_3, int_4, maximum]
self.colour_ramp['colours'] = colours
self.colour_ramp['ranges'] = value_list
area_per_colour = self.get_area(raster_path)
percentage_0 = round(((area_per_colour[colours[0]] / total_area_raster) * 100), 2)
percentage_1 = round(((area_per_colour[colours[1]] / total_area_raster) * 100), 2)
percentage_2 = round(((area_per_colour[colours[2]] / total_area_raster) * 100), 2)
percentage_3 = round(((area_per_colour[colours[3]] / total_area_raster) * 100), 2)
percentage_4 = round(((area_per_colour[colours[4]] / total_area_raster) * 100), 2)
percentage_5 = round(((area_per_colour[colours[5]] / total_area_raster) * 100), 2)
self.colour_ramp['ramp'] = (
QgsColorRampShader.ColorRampItem(value_list[0], QColor(colours[0]), lbl=str('%d (%s ha - %s %%)' % (
value_list[0], area_per_colour[colours[0]], percentage_0))),
QgsColorRampShader.ColorRampItem(value_list[1], QColor(colours[1]),
lbl='%d - %d (%s ha - %s %%)' % (
value_list[0] + 1,
value_list[1],
area_per_colour[colours[1]], percentage_1)),
QgsColorRampShader.ColorRampItem(value_list[2], QColor(colours[2]),
lbl='%d - %d (%s ha - %s %%)' % (
value_list[1] + 1,
value_list[2],
area_per_colour[colours[2]], percentage_2)),
QgsColorRampShader.ColorRampItem(value_list[3], QColor(colours[3]),
lbl='%d - %d (%s ha - %s %%)' % (
value_list[2] + 1,
value_list[3],
area_per_colour[colours[3]], percentage_3)),
QgsColorRampShader.ColorRampItem(value_list[4], QColor(colours[4]),
lbl='%d - %d (%s ha - %s %%)' % (
value_list[3] + 1,
value_list[4],
area_per_colour[colours[4]], percentage_4)),
QgsColorRampShader.ColorRampItem(value_list[5], QColor(colours[5]),
lbl='%d - %d (%s ha - %s %%)' % (
value_list[4] + 1,
value_list[5],
area_per_colour[colours[5]], percentage_5)))
self.mini_style(layer)
# function to determine utm zones
def crs_zone(self, raster_path):
data = gdal.Open(raster_path, gdal.GA_ReadOnly)
if data is None:
print("Raster layer not valid")
else:
projection = osr.SpatialReference(wkt=data.GetProjection())
epsg_code = projection.GetAttrValue('AUTHORITY', 1)
return epsg_code
# Main function that has all the logic
def run(self):
raster = '/tmp/S34E020.hgt'
eps_code = self.crs_zone(raster)
if eps_code == '4326':
projected_raster = processing.run("gdal:warpreproject", {'INPUT': '%s' % raster, 'SOURCE_CRS': None,
'TARGET_CRS': QgsCoordinateReferenceSystem(
'EPSG:3857'), 'RESAMPLING': 0, 'NODATA': None,
'TARGET_RESOLUTION': None,
'OPTIONS': 'COMPRESS=DEFLATE|PREDICTOR=2|ZLEVEL=9',
'DATA_TYPE': 0, 'TARGET_EXTENT': None,
'TARGET_EXTENT_CRS': None, 'MULTITHREADING': True,
'EXTRA': '',
'OUTPUT': 'TEMPORARY_OUTPUT'})
raster_path = projected_raster.get('OUTPUT')
if not raster_path:
raise Exception('Raster was not reprojected')
else:
raster_path = raster
print("The raster layer that is generated is:" + raster_path)
self.style_raster(raster_path)
self.get_area(raster_path)
reporter = ReportClass()
reporter.run()
""
@NyakudyaA
Copy link
Author

should work for both, mind you this was written for QGIS2, you might need to port it to QGIS 3

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment