### Navigation

- index
- modules
- Home »
- Contents »
- Module code »

# 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

    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])])

        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

### Navigation

- index
- modules
- Home »
- Contents »
- Module code »