Source code for safe.impact_functions.volcanic.volcano_building_impact

import numpy
from safe.impact_functions.core import (FunctionProvider,
                                        get_hazard_layer,
                                        get_exposure_layer,
                                        get_question,
                                        format_int)
from safe.storage.vector import Vector
from safe.common.utilities import ugettext as tr
from safe.common.tables import Table, TableRow
from safe.engine.interpolation import (assign_hazard_values_to_exposure_data,
                                       make_circular_polygon)
from safe.common.exceptions import InaSAFEError
from third_party.odict import OrderedDict


[docs]class VolcanoBuildingImpact(FunctionProvider): """Risk plugin for volvano evacuation :author AIFDR :rating 4 :param requires category=='hazard' and \ subcategory in ['volcano'] and \ layertype=='vector' :param requires category=='exposure' and \ subcategory=='structure' and \ layertype=='vector' """ title = tr('Be affected') target_field = 'buildings' parameters = OrderedDict([ ('distances [km]', [3, 5, 10])])
[docs] def run(self, layers): """Risk plugin for flood population evacuation Input layers: List of layers expected to contain H: Raster layer of volcano depth P: Raster layer of population data on the same grid as H Counts number of people exposed to each volcano hazard zones. Return Map of population exposed to volcanic hazard zones Table with number of buildings affected """ # Identify hazard and exposure layers H = get_hazard_layer(layers) # Flood inundation E = get_exposure_layer(layers) is_point_data = False question = get_question(H.get_name(), E.get_name(), self) # Input checks if not H.is_vector: msg = ('Input hazard %s was not a vector layer as expected ' % H.get_name()) raise Exception(msg) msg = ('Input hazard must be a polygon or point layer. ' 'I got %s with layer ' 'type %s' % (H.get_name(), H.get_geometry_name())) if not (H.is_polygon_data or H.is_point_data): raise Exception(msg) if H.is_point_data: # Use concentric circles radii = self.parameters['distances [km]'] is_point_data = True centers = H.get_geometry() attributes = H.get_data() rad_m = [x * 1000 for x in radii] # Convert to meters Z = make_circular_polygon(centers, rad_m, attributes=attributes) # NOTE (Sunni) : I commented out this one because there will be # a permission problem on windows # Z.write_to_file('Marapi_evac_zone_%s.shp' % str(rad_m)) # To check category_title = 'Radius' H = Z #category_names = ['%s m' % x for x in radii] category_names = rad_m name_attribute = 'NAME' # As in e.g. the Smithsonian dataset else: # Use hazard map category_title = 'KRB' # FIXME (Ole): Change to English and use translation system category_names = ['Kawasan Rawan Bencana III', 'Kawasan Rawan Bencana II', 'Kawasan Rawan Bencana I'] name_attribute = 'GUNUNG' # As in e.g. BNPB hazard map # Get names of volcanos considered if name_attribute in H.get_attribute_names(): D = {} for att in H.get_data(): # Run through all polygons and get unique names D[att[name_attribute]] = None volcano_names = '' for name in D: volcano_names += '%s, ' % name volcano_names = volcano_names[:-2] # Strip trailing ', ' else: volcano_names = tr('Not specified in data') if not category_title in H.get_attribute_names(): msg = ('Hazard data %s did not contain expected ' 'attribute %s ' % (H.get_name(), category_title)) raise InaSAFEError(msg) # Run interpolation function for polygon2raster P = assign_hazard_values_to_exposure_data(H, E) # Initialise attributes of output dataset with all attributes # from input polygon and a population count of zero new_attributes = H.get_data() categories = {} for attr in new_attributes: attr[self.target_field] = 0 cat = attr[category_title] categories[cat] = 0 # Count affected population per polygon and total # total_affected = 0 for attr in P.get_data(): # Update building count for associated polygon poly_id = attr['polygon_id'] if poly_id is not None: new_attributes[poly_id][self.target_field] += 1 # Update building count for each category cat = new_attributes[poly_id][category_title] categories[cat] += 1 # Update total # total_affected += 1 # Count totals total = len(E) # Generate simple impact report blank_cell = '' table_body = [question, TableRow([tr('Volcanos considered'), '%s' % volcano_names, blank_cell], header=True), TableRow([tr('Distance [km]'), tr('Total'), tr('Cumulative')], header=True)] cum = 0 for name in category_names: # prevent key error count = categories.get(name, 0) cum += count if is_point_data: name = int(name) / 1000 table_body.append(TableRow([name, format_int(count), format_int(cum)])) table_body.append(TableRow(tr('Map shows buildings affected in ' 'each of volcano hazard polygons.'))) impact_table = Table(table_body).toNewlineFreeString() # Extend impact report for on-screen display table_body.extend([TableRow(tr('Notes'), header=True), tr('Total number of buildings %s in the viewable ' 'area') % format_int(total), tr('Only buildings available in OpenStreetMap ' 'are considered.')]) impact_summary = Table(table_body).toNewlineFreeString() map_title = tr('Buildings affected by volcanic hazard zone') # Define classes for legend for flooded building counts colours = ['#FFFFFF', '#38A800', '#79C900', '#CEED00', '#FFCC00', '#FF6600', '#FF0000', '#7A0000'] building_counts = [x[self.target_field] for x in new_attributes] cls = [0] + numpy.linspace(1, max(building_counts), len(colours)).tolist() # Define style info for output polygons showing population counts style_classes = [] for i, colour in enumerate(colours): lo = cls[i] hi = cls[i + 1] if i == 0: label = tr('0') else: label = tr('%i - %i') % (lo, hi) entry = dict(label=label, colour=colour, min=lo, max=hi, transparency=0, size=1) style_classes.append(entry) # Override style info with new classes and name style_info = dict(target_field=self.target_field, style_classes=style_classes, legend_title=tr('Building Count')) # Create vector layer and return V = Vector(data=new_attributes, projection=H.get_projection(), geometry=H.get_geometry(as_geometry_objects=True), name=tr('Buildings affected by volcanic hazard zone'), keywords={'impact_summary': impact_summary, 'impact_table': impact_table, 'map_title': map_title, 'target_field': self.target_field}, style_info=style_info) return V