Last active
March 24, 2023 04:29
-
-
Save NyakudyaA/b4640ec9d2b5f43fa456083b61cfd12f to your computer and use it in GitHub Desktop.
PyQGIS Implementation of generating raster statistics
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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() | |
| "" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
should work for both, mind you this was written for QGIS2, you might need to port it to QGIS 3