Python interpolation.assign_hazard_values_to_exposure_data函数代码示例

    def run(layers):
        """Risk plugin for earthquake school damage
        """

        # Extract data
        H = get_hazard_layer(layers)    # Ground shaking
        E = get_exposure_layer(layers)  # Building locations

        # Interpolate hazard level to building locations
        H = assign_hazard_values_to_exposure_data(H, E,
                                             attribute_name='MMI')

        # Extract relevant numerical data
        coordinates = E.get_geometry()
        shaking = H.get_data()

        # Calculate building damage
        building_damage = []
        for i in range(len(shaking)):
            x = float(shaking[i]['MMI'])
            if x < 6.0 or (x != x):  # x != x -> check for nan pre python 2.6
                value = 0.0
            else:
                value = (0.692 * (x ** 4) -
                         15.82 * (x ** 3) +
                         135.0 * (x ** 2) -
                         509.0 * x + 714.4)

            building_damage.append({'DAMAGE': value, 'MMI': x})

        # Create new layer and return
        V = Vector(data=building_damage,
                   projection=E.get_projection(),
                   geometry=coordinates)
        return V

    def run(self, layers=None):
        """Risk plugin for classified polygon hazard on building/structure.

        Counts number of building exposed to each hazard zones.

        :param layers: List of layers expected to contain.
                * hazard_layer: Hazard layer
                * exposure_layer: Vector layer of structure data on
                the same grid as hazard_layer

        :returns: Map of building exposed to each hazard zones.
                  Table with number of buildings affected
        :rtype: dict
        """
        self.validate()
        self.prepare(layers)

        # Target Field
        target_field = 'zone'

        # Not affected string in the target field
        not_affected_value = 'Not Affected'

        # Parameters
        hazard_zone_attribute = self.parameters['hazard zone attribute']

        # Identify hazard and exposure layers
        hazard_layer = self.hazard
        exposure_layer = self.exposure

        # Input checks
        if not hazard_layer.is_polygon_data:
            message = (
                'Input hazard must be a polygon. I got %s with '
                'layer type %s' %
                (hazard_layer.get_name(), hazard_layer.get_geometry_name()))
            raise Exception(message)

        # Check if hazard_zone_attribute exists in hazard_layer
        if hazard_zone_attribute not in hazard_layer.get_attribute_names():
            message = (
                'Hazard data %s does not contain expected attribute %s ' %
                (hazard_layer.get_name(), hazard_zone_attribute))
            # noinspection PyExceptionInherit
            raise InaSAFEError(message)

        # Find the target field name that has no conflict with default
        # target
        attribute_names = hazard_layer.get_attribute_names()
        target_field = get_non_conflicting_attribute_name(
            target_field, attribute_names)

        # Hazard zone categories from hazard layer
        self.hazard_zones = list(
            set(hazard_layer.get_data(hazard_zone_attribute)))

        self.buildings = {}
        self.affected_buildings = OrderedDict()
        for hazard_zone in self.hazard_zones:
            self.affected_buildings[hazard_zone] = {}

        # Run interpolation function for polygon2polygon
        interpolated_layer = assign_hazard_values_to_exposure_data(
            hazard_layer, exposure_layer, attribute_name=None)

        # Extract relevant interpolated data
        attribute_names = interpolated_layer.get_attribute_names()
        features = interpolated_layer.get_data()

        for i in range(len(features)):
            hazard_value = features[i][hazard_zone_attribute]
            if not hazard_value:
                hazard_value = not_affected_value
            features[i][target_field] = hazard_value
            usage = get_osm_building_usage(attribute_names, features[i])
            if usage is None:
                usage = tr('Unknown')
            if usage not in self.buildings:
                self.buildings[usage] = 0
                for category in self.affected_buildings.keys():
                    self.affected_buildings[category][usage] = OrderedDict(
                        [(tr('Buildings Affected'), 0)])
            self.buildings[usage] += 1
            if hazard_value in self.affected_buildings.keys():
                self.affected_buildings[hazard_value][usage][
                    tr('Buildings Affected')] += 1

        # Lump small entries and 'unknown' into 'other' category
        self._consolidate_to_other()

        # Generate simple impact report
        impact_summary = impact_table = self.generate_html_report()

        # Create style
        categories = self.hazard_zones
        categories.append(not_affected_value)
        colours = color_ramp(len(categories))
        style_classes = []

        i = 0
#.........这里部分代码省略.........

    def run(self, layers):
        """Risk plugin for volcano hazard on building/structure

        Input
          layers: List of layers expected to contain
              my_hazard: Hazard layer of volcano
              my_exposure: Vector layer of structure data on
              the same grid as my_hazard

        Counts number of building exposed to each volcano hazard zones.

        Return
          Map of building exposed to volcanic hazard zones
          Table with number of buildings affected
        """

        # Identify hazard and exposure layers
        my_hazard = get_hazard_layer(layers)  # Volcano hazard layer
        my_exposure = get_exposure_layer(layers)
        is_point_data = False

        question = get_question(my_hazard.get_name(), my_exposure.get_name(), self)

        # Input checks
        if not my_hazard.is_vector:
            msg = "Input hazard %s  was not a vector layer as expected " % my_hazard.get_name()
            raise Exception(msg)

        msg = "Input hazard must be a polygon or point layer. I got %s " "with layer type %s" % (
            my_hazard.get_name(),
            my_hazard.get_geometry_name(),
        )
        if not (my_hazard.is_polygon_data or my_hazard.is_point_data):
            raise Exception(msg)

        if my_hazard.is_point_data:
            # Use concentric circles
            radii = self.parameters["distances [km]"]
            is_point_data = True

            centers = my_hazard.get_geometry()
            attributes = my_hazard.get_data()
            rad_m = [x * 1000 for x in radii]  # Convert to meters
            Z = make_circular_polygon(centers, rad_m, attributes=attributes)
            # To check
            category_title = "Radius"
            my_hazard = Z

            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 my_hazard.get_attribute_names():
            D = {}
            for att in my_hazard.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 my_hazard.get_attribute_names():
            msg = "Hazard data %s did not contain expected " "attribute %s " % (my_hazard.get_name(), category_title)
            # noinspection PyExceptionInherit
            raise InaSAFEError(msg)

        # Run interpolation function for polygon2raster
        P = assign_hazard_values_to_exposure_data(my_hazard, my_exposure)

        # Initialise attributes of output dataset with all attributes
        # from input polygon and a building count of zero
        new_attributes = my_hazard.get_data()

        categories = {}
        for attr in new_attributes:
            attr[self.target_field] = 0
            cat = attr[category_title]
            categories[cat] = 0

        # Count impacted building per polygon and total
        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]
#.........这里部分代码省略.........

    def run(self, layers):
        """Risk plugin for Padang building survey
        """

        # Extract data
        H = get_hazard_layer(layers)  # Ground shaking
        E = get_exposure_layer(layers)  # Building locations

        datatype = E.get_keywords()["datatype"]
        vclass_tag = "ITB_Class"
        if datatype.lower() == "osm":
            # Map from OSM attributes to the ITB building classes
            #            Emap = osm2itb(E)
            print "osm2itb has not been implemented"
        elif datatype.lower() == "sigab":
            #            Emap = sigabitb(E)
            print "sigab2itb has not been implemented"
        elif datatype.lower() == "itb":
            Emap = E

        # Interpolate hazard level to building locations
        Hi = assign_hazard_values_to_exposure_data(H, Emap, attribute_name="MMI")

        # Extract relevant numerical data
        coordinates = Emap.get_geometry()
        shaking = Hi.get_data()
        N = len(shaking)

        # List attributes to carry forward to result layer
        attributes = Emap.get_attribute_names()
        # Calculate building damage
        count50 = 0
        count25 = 0
        count10 = 0
        count0 = 0
        building_damage = []
        for i in range(N):
            mmi = float(shaking[i]["MMI"])

            building_class = Emap.get_data(vclass_tag, i)

            building_type = str(building_class)
            damage_params = vul_curves[building_type]
            beta = damage_params["beta"]
            median = damage_params["median"]

            msg = "Invalid parameter value for " + building_type
            verify(beta + median > 0.0, msg)
            percent_damage = lognormal_cdf(mmi, median=median, sigma=beta) * 100

            # Collect shake level and calculated damage
            result_dict = {self.target_field: percent_damage, "MMI": mmi}

            # Carry all orginal attributes forward
            for key in attributes:
                result_dict[key] = Emap.get_data(key, i)

            # Record result for this feature
            building_damage.append(result_dict)

            # Debugging
            # if percent_damage > 0.01:
            #    print mmi, percent_damage

            # Calculate statistics
            if percent_damage < 10:
                count0 += 1

            if 10 <= percent_damage < 33:
                count10 += 1

            if 33 <= percent_damage < 66:
                count25 += 1

            if 66 <= percent_damage:
                count50 += 1

        #        fid.close()
        # Create report
        Hname = H.get_name()
        Ename = E.get_name()
        impact_summary = '<b>In case of "%s" the estimated impact to ' '"%s" ' "is&#58;</b><br><br><p>" % (Hname, Ename)
        impact_summary += (
            '<table border="0" width="320px">'
            "   <tr><th><b>%s</b></th><th><b>%s</b></th></th>"
            "   <tr></tr>"
            "   <tr><td>%s&#58;</td><td>%i</td></tr>"
            "   <tr><td>%s (<10%%)&#58;</td><td>%i</td></tr>"
            "   <tr><td>%s (10-33%%)&#58;</td><td>%i</td></tr>"
            "   <tr><td>%s (33-66%%)&#58;</td><td>%i</td></tr>"
            "   <tr><td>%s (66-100%%)&#58;</td><td>%i</td></tr>"
            "</table></font>"
            % (
                tr("Buildings"),
                tr("Total"),
                tr("All"),
                N,
                tr("No damage"),
                count0,
                tr("Low damage"),
#.........这里部分代码省略.........

    def run(self):
        """Run volcano point population evacuation Impact Function.

        Counts number of people exposed to volcano event.

        :returns: Map of population exposed to the volcano hazard zone.
            The returned dict will include a table with number of people
            evacuated and supplies required.
        :rtype: dict

        :raises:
            * Exception - When hazard layer is not vector layer
            * RadiiException - When radii are not valid (they need to be
                monotonically increasing)
        """

        # Parameters
        radii = self.parameters['distances'].value

        # Get parameters from layer's keywords
        volcano_name_attribute = self.hazard.keyword('volcano_name_field')

        data_table = self.hazard.layer.get_data()

        # Get names of volcanoes considered
        if volcano_name_attribute in self.hazard.layer.get_attribute_names():
            volcano_name_list = []
            # Run through all polygons and get unique names
            for row in data_table:
                volcano_name_list.append(row[volcano_name_attribute])

            volcano_names = ''
            for radius in volcano_name_list:
                volcano_names += '%s, ' % radius
            self.volcano_names = volcano_names[:-2]  # Strip trailing ', '

        # Run interpolation function for polygon2raster
        interpolated_layer, covered_exposure_layer = \
            assign_hazard_values_to_exposure_data(
                self.hazard.layer,
                self.exposure.layer,
                attribute_name=self.target_field
            )

        # Initialise affected population per categories
        for radius in radii:
            category = 'Radius %s km ' % format_int(radius)
            self.affected_population[category] = 0

        if has_no_data(self.exposure.layer.get_data(nan=True)):
            self.no_data_warning = True
        # Count affected population per polygon and total
        for row in interpolated_layer.get_data():
            # Get population at this location
            population = row[self.target_field]
            if not numpy.isnan(population):
                population = float(population)
                # Update population count for this category
                category = 'Radius %s km ' % format_int(
                    row[self.hazard_zone_attribute])
                self.affected_population[category] += population

        # Count totals
        self.total_population = population_rounding(
            int(numpy.nansum(self.exposure.layer.get_data())))

        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]

        # Create style
        colours = ['#FFFFFF', '#38A800', '#79C900', '#CEED00',
                   '#FFCC00', '#FF6600', '#FF0000', '#7A0000']
        classes = create_classes(
            covered_exposure_layer.get_data().flat[:], len(colours))
        interval_classes = humanize_class(classes)
        # Define style info for output polygons showing population counts
        style_classes = []
        for i in xrange(len(colours)):
            style_class = dict()
            style_class['label'] = create_label(interval_classes[i])
            if i == 1:
                label = create_label(
                    interval_classes[i],
                    tr('Low Population [%i people/cell]' % classes[i]))
            elif i == 4:
                label = create_label(
                    interval_classes[i],
                    tr('Medium Population [%i people/cell]' % classes[i]))
            elif i == 7:
                label = create_label(
                    interval_classes[i],
                    tr('High Population [%i people/cell]' % classes[i]))
            else:
                label = create_label(interval_classes[i])

            style_class['label'] = label
            style_class['quantity'] = classes[i]
            style_class['colour'] = colours[i]
#.........这里部分代码省略.........

    def run(self, layers):
        """Risk plugin for earthquake school damage
        """

        # Extract data
        H = get_hazard_layer(layers)    # Ground shaking
        E = get_exposure_layer(layers)  # Building locations

        keywords = E.get_keywords()
        if 'datatype' in keywords:
            datatype = keywords['datatype']
            if datatype.lower() == 'osm':
                # Map from OSM attributes to the guideline classes (URM and RM)
                E = osm2bnpb(E, target_attribute=self.vclass_tag)
            elif datatype.lower() == 'sigab':
                # Map from SIGAB attributes to the guideline classes
                # (URM and RM)
                E = sigab2bnpb(E)
            else:
                E = unspecific2bnpb(E, target_attribute=self.vclass_tag)
        else:
            E = unspecific2bnpb(E, target_attribute=self.vclass_tag)

        # Interpolate hazard level to building locations
        H = assign_hazard_values_to_exposure_data(H, E,
                                             attribute_name='MMI')

        # Extract relevant numerical data
        coordinates = E.get_geometry()
        shaking = H.get_data()
        N = len(shaking)

        # List attributes to carry forward to result layer
        attributes = E.get_attribute_names()

        # Calculate building damage
        count3 = 0
        count2 = 0
        count1 = 0
        count_unknown = 0
        building_damage = []
        for i in range(N):
            mmi = float(shaking[i]['MMI'])

            building_class = E.get_data(self.vclass_tag, i)
            lo, hi = damage_parameters[building_class]

            if numpy.isnan(mmi):
                # If we don't know the shake level assign Not-a-Number
                damage = numpy.nan
                count_unknown += 1
            elif mmi < lo:
                damage = 1  # Low
                count1 += 1
            elif lo <= mmi < hi:
                damage = 2  # Medium
                count2 += 1
            elif mmi >= hi:
                damage = 3  # High
                count3 += 1
            else:
                msg = 'Undefined shakelevel %s' % str(mmi)
                raise Exception(msg)

            # Collect shake level and calculated damage
            result_dict = {self.target_field: damage,
                           'MMI': mmi}

            # Carry all orginal attributes forward
            for key in attributes:
                result_dict[key] = E.get_data(key, i)

            # Record result for this feature
            building_damage.append(result_dict)

        # Create report
        impact_summary = ('<table border="0" width="320px">'
                   '   <tr><th><b>%s</b></th><th><b>%s</b></th></th>'
                    '   <tr></tr>'
                    '   <tr><td>%s&#58;</td><td>%s</td></tr>'
                    '   <tr><td>%s (10-25%%)&#58;</td><td>%s</td></tr>'
                    '   <tr><td>%s (25-50%%)&#58;</td><td>%s</td></tr>'
                    '   <tr><td>%s (50-100%%)&#58;</td><td>%s</td></tr>'
                    % (tr('Buildings'), tr('Total'),
                       tr('All'), format_int(N),
                       tr('Low damage'), format_int(count1),
                       tr('Medium damage'), format_int(count2),
                       tr('High damage'), format_int(count3)))
        impact_summary += ('   <tr><td>%s (NaN)&#58;</td><td>%s</td></tr>'
                    % ('Unknown', format_int(count_unknown)))
        impact_summary += '</table>'

        # Create style
        style_classes = [dict(label=tr('Low damage'), min=0.5, max=1.5,
                              colour='#fecc5c', transparency=0),
                         dict(label=tr('Medium damage'), min=1.5, max=2.5,
                              colour='#fd8d3c', transparency=0),
                         dict(label=tr('High damage'), min=2.5, max=3.5,
                              colour='#f31a1c', transparency=0)]
        style_info = dict(target_field=self.target_field,
#.........这里部分代码省略.........

    def run(self, layers):
        """Flood impact to buildings (e.g. from Open Street Map)
        """

        # Extract data
        H = get_hazard_layer(layers)    # Depth
        E = get_exposure_layer(layers)  # Building locations

        question = get_question(H.get_name(),
                                E.get_name(),
                                self)

        # Interpolate hazard level to building locations
        I = assign_hazard_values_to_exposure_data(H, E)

        # Extract relevant exposure data
        #attribute_names = I.get_attribute_names()
        attributes = I.get_data()
        N = len(I)

        # Calculate road impact
        count = 0
        #flooded_distance = 0
        for i in range(N):
            # Use interpolated polygon attribute
            atts = attributes[i]

            if 'FLOODPRONE' in atts:
                res = atts['FLOODPRONE']
                if res is None:
                    x = False
                else:
                    x = res.lower() == 'yes'
            else:
                # If there isn't a flood prone attribute,
                # assume that building is wet if inside polygon
                # as flag by generic attribute AFFECTED
                res = atts['Affected']
                if res is None:
                    x = False
                else:
                    x = res

            # Count all roads
            if x is True:
                # Count total affected roads
                count += 1

            # Add calculated impact to existing attributes
            attributes[i][self.target_field] = x
            if i == 0:
                print attributes[0].keys()

        # Generate simple impact report
        table_body = [question,
                      TableRow([tr('Building type'),
                                tr('Temporarily closed'),
                                tr('Total')],
                               header=True),
                      TableRow([tr('All'), count, N])]
        impact_summary = Table(table_body).toNewlineFreeString()
        #impact_table = impact_summary
        map_title = tr('Roads inundated')

        # Create style
        style_classes = [dict(label=tr('Not Flooded'), min=0, max=0,
                              colour='#1EFC7C', transparency=0, size=1),
                         dict(label=tr('Flooded'), min=1, max=1,
                              colour='#F31A1C', transparency=0, size=1)]
        style_info = dict(target_field=self.target_field,
                          style_classes=style_classes)

        # Create vector layer and return
        V = Vector(data=attributes,
                   projection=I.get_projection(),
                   geometry=I.get_geometry(),
                   geometry_type=I.get_geometry_type(),
                   name=tr('Estimated roads affected'),
                   keywords={'impact_summary': impact_summary,
                             'map_title': map_title,
                             'target_field': self.target_field},
                   style_info=style_info)
        return V

    def run(self, layers):
        """Risk plugin for volcano population evacuation

        Input
          layers: List of layers expected to contain
              H: Vector polygon layer of volcano impact zones
              P: Raster layer of population data on the same grid as H

        Counts number of people exposed to volcano event.

        Return
          Map of population exposed to the volcano hazard zone.
          Table with number of people evacuated and supplies required.
        """

        # Identify hazard and exposure layers
        H = get_hazard_layer(layers)  # Flood inundation
        E = get_exposure_layer(layers)

        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['distance [km]']

            centers = H.get_geometry()
            attributes = H.get_data()
            rad_m = [x * 1000 for x in radii]  # Convert to meters
            H = make_circular_polygon(centers,
                                      rad_m,
                                      attributes=attributes)
            # NOTE (Sunni) : I commented out this one because there will be
            # a permission problem on windows
            #H.write_to_file('Evac_zones_%s.shp' % str(radii))  # To check

            category_title = 'Radius'
            category_header = tr('Distance [km]')
            category_names = radii

            name_attribute = 'NAME'  # As in e.g. the Smithsonian dataset
        else:
            # Use hazard map
            category_title = 'KRB'
            category_header = tr('Category')

            # 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
            attributes = H.get_data()

        # 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,
                                                  attribute_name='population')

        # 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
#.........这里部分代码省略.........

    def run(self):
        """Risk plugin for flood population evacuation.

        Counts number of people exposed to areas identified as flood prone

        :returns: Map of population exposed to flooding Table with number of
            people evacuated and supplies required.
        :rtype: tuple
        """
        self.validate()
        self.prepare()

        self.provenance.append_step(
            'Calculating Step',
            'Impact function is calculating the impact.')

        # Get parameters from layer's keywords
        self.hazard_class_attribute = self.hazard.keyword('field')
        self.hazard_class_mapping = self.hazard.keyword('value_map')

        # Get the IF parameters
        self._evacuation_percentage = (
            self.parameters['evacuation_percentage'].value)

        # Check that hazard is polygon type
        if not self.hazard.layer.is_polygon_data:
            message = (
                'Input hazard must be a polygon layer. I got %s with layer '
                'type %s' % (
                    self.hazard.name,
                    self.hazard.layer.get_geometry_name()))
            raise Exception(message)

        if has_no_data(self.exposure.layer.get_data(nan=True)):
            self.no_data_warning = True

        # Check that affected field exists in hazard layer
        if (self.hazard_class_attribute in
                self.hazard.layer.get_attribute_names()):
            self.use_affected_field = True

        # Run interpolation function for polygon2raster
        interpolated_layer, covered_exposure = \
            assign_hazard_values_to_exposure_data(
                self.hazard.layer,
                self.exposure.layer,
                attribute_name=self.target_field)

        # Data for manipulating the covered_exposure layer
        new_covered_exposure_data = covered_exposure.get_data()
        covered_exposure_top_left = numpy.array([
            covered_exposure.get_geotransform()[0],
            covered_exposure.get_geotransform()[3]])
        covered_exposure_dimension = numpy.array([
            covered_exposure.get_geotransform()[1],
            covered_exposure.get_geotransform()[5]])

        # Count affected population per polygon, per category and total
        total_affected_population = 0
        for attr in interpolated_layer.get_data():
            affected = False
            if self.use_affected_field:
                row_affected_value = attr[self.hazard_class_attribute]
                if row_affected_value is not None:
                    affected = get_key_for_value(
                        row_affected_value, self.hazard_class_mapping)
            else:
                # assume that every polygon is affected (see #816)
                affected = self.wet

            if affected == self.wet:
                # Get population at this location
                population = attr[self.target_field]
                if not numpy.isnan(population):
                    population = float(population)
                    total_affected_population += population
            else:
                # If it's not affected, set the value of the impact layer to 0
                grid_point = attr['grid_point']
                index = numpy.floor(
                    (grid_point - covered_exposure_top_left) / (
                        covered_exposure_dimension)).astype(int)
                new_covered_exposure_data[index[1]][index[0]] = 0

        # Estimate number of people in need of evacuation
        if self.use_affected_field:
            affected_population = tr(
                'People within hazard field ("%s") of value "%s"') % (
                    self.hazard_class_attribute,
                    ','.join([
                        unicode(hazard_class) for
                        hazard_class in self.hazard_class_mapping[self.wet]
                    ]))
        else:
            affected_population = tr('People within any hazard polygon.')

        self.affected_population[affected_population] = (
            total_affected_population)

        self.total_population = int(
#.........这里部分代码省略.........

    def run(self, layers):
        """Risk plugin for volcano population evacuation.

        :param layers: List of layers expected to contain where two layers
            should be present.

            * hazard_layer: Vector polygon layer of volcano impact zones
            * exposure_layer: Raster layer of population data on the same grid
                as hazard_layer

        Counts number of people exposed to volcano event.

        :returns: Map of population exposed to the volcano hazard zone.
            The returned dict will include a table with number of people
            evacuated and supplies required.
        :rtype: dict

        :raises:
            * Exception - When hazard layer is not vector layer
            * RadiiException - When radii are not valid (they need to be
                monotonically increasing)
        """

        # Identify hazard and exposure layers
        hazard_layer = get_hazard_layer(layers)  # Volcano KRB
        exposure_layer = get_exposure_layer(layers)

        question = get_question(
            hazard_layer.get_name(), exposure_layer.get_name(), self)

        # Input checks
        if not hazard_layer.is_vector:
            msg = ('Input hazard %s  was not a vector layer as expected '
                   % hazard_layer.get_name())
            raise Exception(msg)

        msg = ('Input hazard must be a polygon or point layer. I got %s with '
               'layer type %s' % (hazard_layer.get_name(),
                                  hazard_layer.get_geometry_name()))
        if not (hazard_layer.is_polygon_data or hazard_layer.is_point_data):
            raise Exception(msg)

        data_table = hazard_layer.get_data()
        if hazard_layer.is_point_data:
            # Use concentric circles
            radii = self.parameters['distance [km]']

            centers = hazard_layer.get_geometry()
            rad_m = [x * 1000 for x in radii]  # Convert to meters
            hazard_layer = buffer_points(centers, rad_m, data_table=data_table)

            category_title = 'Radius'
            category_header = tr('Distance [km]')
            category_names = radii

            name_attribute = 'NAME'  # As in e.g. the Smithsonian dataset
        else:
            # Use hazard map
            category_title = 'KRB'
            category_header = tr('Category')

            # 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 volcanoes considered
        if name_attribute in hazard_layer.get_attribute_names():
            volcano_name_list = []
            # Run through all polygons and get unique names
            for row in data_table:
                volcano_name_list.append(row[name_attribute])

            volcano_names = ''
            for name in volcano_name_list:
                volcano_names += '%s, ' % name
            volcano_names = volcano_names[:-2]  # Strip trailing ', '
        else:
            volcano_names = tr('Not specified in data')

        # Check if category_title exists in hazard_layer
        if not category_title in hazard_layer.get_attribute_names():
            msg = ('Hazard data %s did not contain expected '
                   'attribute %s ' % (hazard_layer.get_name(), category_title))
            # noinspection PyExceptionInherit
            raise InaSAFEError(msg)

        # Find the target field name that has no conflict with default target
        attribute_names = hazard_layer.get_attribute_names()
        new_target_field = get_non_conflicting_attribute_name(
            self.target_field, attribute_names)
        self.target_field = new_target_field

        # Run interpolation function for polygon2raster
        interpolated_layer = assign_hazard_values_to_exposure_data(
            hazard_layer, exposure_layer, attribute_name=self.target_field)

        # Initialise data_table of output dataset with all data_table
#.........这里部分代码省略.........

    def run(self, layers):
        """Flood impact to buildings (e.g. from Open Street Map).

         :param layers: List of layers expected to contain.
                * hazard_layer: Hazard layer of flood
                * exposure_layer: Vector layer of structure data on
                the same grid as hazard_layer
        """
        threshold = self.parameters['threshold [m]']  # Flood threshold [m]

        verify(isinstance(threshold, float),
               'Expected thresholds to be a float. Got %s' % str(threshold))

        # Extract data
        hazard_layer = get_hazard_layer(layers)  # Depth
        exposure_layer = get_exposure_layer(layers)  # Building locations

        question = get_question(
            hazard_layer.get_name(),
            exposure_layer.get_name(),
            self)

        # Determine attribute name for hazard levels
        if hazard_layer.is_raster:
            mode = 'grid'
            hazard_attribute = 'depth'
        else:
            mode = 'regions'
            hazard_attribute = None

        # Interpolate hazard level to building locations
        interpolated_layer = assign_hazard_values_to_exposure_data(
            hazard_layer, exposure_layer, attribute_name=hazard_attribute)

        # Extract relevant exposure data
        attribute_names = interpolated_layer.get_attribute_names()
        features = interpolated_layer.get_data()
        total_features = len(interpolated_layer)
        buildings = {}

        # The number of affected buildings
        affected_count = 0

        # The variable for grid mode
        inundated_count = 0
        wet_count = 0
        dry_count = 0
        inundated_buildings = {}
        wet_buildings = {}
        dry_buildings = {}

        # The variable for regions mode
        affected_buildings = {}

        if mode == 'grid':
            for i in range(total_features):
                # Get the interpolated depth
                water_depth = float(features[i]['depth'])
                if water_depth <= 0:
                    inundated_status = 0  # dry
                elif water_depth >= threshold:
                    inundated_status = 1  # inundated
                else:
                    inundated_status = 2  # wet

                # Count affected buildings by usage type if available
                usage = get_osm_building_usage(attribute_names, features[i])
                if usage is not None and usage != 0:
                    key = usage
                else:
                    key = 'unknown'

                if key not in buildings:
                    buildings[key] = 0
                    inundated_buildings[key] = 0
                    wet_buildings[key] = 0
                    dry_buildings[key] = 0

                # Count all buildings by type
                buildings[key] += 1
                if inundated_status is 0:
                    # Count dry buildings by type
                    dry_buildings[key] += 1
                    # Count total dry buildings
                    dry_count += 1
                if inundated_status is 1:
                    # Count inundated buildings by type
                    inundated_buildings[key] += 1
                    # Count total dry buildings
                    inundated_count += 1
                if inundated_status is 2:
                    # Count wet buildings by type
                    wet_buildings[key] += 1
                    # Count total wet buildings
                    wet_count += 1
                # Add calculated impact to existing attributes
                features[i][self.target_field] = inundated_status
        elif mode == 'regions':
            for i in range(total_features):
                # Use interpolated polygon attribute
#.........这里部分代码省略.........

    def run(self, layers):
        """Earthquake impact to buildings (e.g. from OpenStreetMap).

        :param layers: All the input layers (Hazard Layer and Exposure Layer)
        """

        LOGGER.debug('Running earthquake building impact')

        # merely initialize
        building_value = 0
        contents_value = 0

        # Thresholds for mmi breakdown.
        t0 = self.parameters['low_threshold']
        t1 = self.parameters['medium_threshold']
        t2 = self.parameters['high_threshold']

        # Class Attribute and Label.

        class_1 = {'label': tr('Low'), 'class': 1}
        class_2 = {'label': tr('Medium'), 'class': 2}
        class_3 = {'label': tr('High'), 'class': 3}

        # Extract data
        hazard_layer = get_hazard_layer(layers)  # Depth
        exposure_layer = get_exposure_layer(layers)  # Building locations

        question = get_question(
            hazard_layer.get_name(),
            exposure_layer.get_name(),
            self
        )

        # Define attribute name for hazard levels.
        hazard_attribute = 'mmi'

        # Determine if exposure data have NEXIS attributes.
        attribute_names = exposure_layer.get_attribute_names()
        if ('FLOOR_AREA' in attribute_names and
            'BUILDING_C' in attribute_names and
                'CONTENTS_C' in attribute_names):
            is_nexis = True
        else:
            is_nexis = False

        # Interpolate hazard level to building locations.
        my_interpolate_result = assign_hazard_values_to_exposure_data(
            hazard_layer,
            exposure_layer,
            attribute_name=hazard_attribute
        )

        # Extract relevant exposure data
        #attribute_names = my_interpolate_result.get_attribute_names()
        attributes = my_interpolate_result.get_data()

        interpolate_size = len(my_interpolate_result)

        # Calculate building impact
        lo = 0
        me = 0
        hi = 0
        building_values = {}
        contents_values = {}
        for key in range(4):
            building_values[key] = 0
            contents_values[key] = 0
        for i in range(interpolate_size):
            # Classify building according to shake level
            # and calculate dollar losses

            if is_nexis:
                try:
                    area = float(attributes[i]['FLOOR_AREA'])
                except (ValueError, KeyError):
                    #print 'Got area', attributes[i]['FLOOR_AREA']
                    area = 0.0

                try:
                    building_value_density = float(attributes[i]['BUILDING_C'])
                except (ValueError, KeyError):
                    #print 'Got bld value', attributes[i]['BUILDING_C']
                    building_value_density = 0.0

                try:
                    contents_value_density = float(attributes[i]['CONTENTS_C'])
                except (ValueError, KeyError):
                    #print 'Got cont value', attributes[i]['CONTENTS_C']
                    contents_value_density = 0.0

                building_value = building_value_density * area
                contents_value = contents_value_density * area

            try:
                x = float(attributes[i][hazard_attribute])  # MMI
            except TypeError:
                x = 0.0
            if t0 <= x < t1:
                lo += 1
                cls = 1
#.........这里部分代码省略.........

    def run(self):
        """Run volcano point population evacuation Impact Function.

        Counts number of people exposed to volcano event.

        :returns: Map of population exposed to the volcano hazard zone.
            The returned dict will include a table with number of people
            evacuated and supplies required.
        :rtype: dict

        :raises:
            * Exception - When hazard layer is not vector layer
            * RadiiException - When radii are not valid (they need to be
                monotonically increasing)
        """
        self.validate()
        self.prepare()

        # Parameters
        radii = self.parameters['distances'].value

        # Get parameters from layer's keywords
        volcano_name_attribute = self.hazard.keyword('volcano_name_field')

        # Input checks
        if not self.hazard.layer.is_point_data:
            msg = (
                'Input hazard must be a polygon or point layer. I got %s with '
                'layer type %s' % (
                    self.hazard.name, self.hazard.layer.get_geometry_name()))
            raise Exception(msg)

        data_table = self.hazard.layer.get_data()

        # Use concentric circles
        category_title = 'Radius'

        centers = self.hazard.layer.get_geometry()
        rad_m = [x * 1000 for x in radii]  # Convert to meters
        hazard_layer = buffer_points(
            centers, rad_m, category_title, data_table=data_table)

        # Get names of volcanoes considered
        if volcano_name_attribute in hazard_layer.get_attribute_names():
            volcano_name_list = []
            # Run through all polygons and get unique names
            for row in data_table:
                volcano_name_list.append(row[volcano_name_attribute])

            volcano_names = ''
            for radius in volcano_name_list:
                volcano_names += '%s, ' % radius
            self.volcano_names = volcano_names[:-2]  # Strip trailing ', '

        # Run interpolation function for polygon2raster
        interpolated_layer, covered_exposure_layer = \
            assign_hazard_values_to_exposure_data(
                hazard_layer,
                self.exposure.layer,
                attribute_name=self.target_field
            )

        # Initialise affected population per categories
        for radius in rad_m:
            category = 'Distance %s km ' % format_int(radius)
            self.affected_population[category] = 0

        if has_no_data(self.exposure.layer.get_data(nan=True)):
            self.no_data_warning = True
        # Count affected population per polygon and total
        for row in interpolated_layer.get_data():
            # Get population at this location
            population = row[self.target_field]
            if not numpy.isnan(population):
                population = float(population)
                # Update population count for this category
                category = 'Distance %s km ' % format_int(
                    row[category_title])
                self.affected_population[category] += population

        # Count totals
        self.total_population = population_rounding(
            int(numpy.nansum(self.exposure.layer.get_data())))

        self.minimum_needs = [
            parameter.serialize() for parameter in
            filter_needs_parameters(self.parameters['minimum needs'])
        ]

        impact_table = impact_summary = self.html_report()

        # Create style
        colours = ['#FFFFFF', '#38A800', '#79C900', '#CEED00',
                   '#FFCC00', '#FF6600', '#FF0000', '#7A0000']
        classes = create_classes(
            covered_exposure_layer.get_data().flat[:], len(colours))
        interval_classes = humanize_class(classes)
        # Define style info for output polygons showing population counts
        style_classes = []
        for i in xrange(len(colours)):
#.........这里部分代码省略.........

    def run(self, layers):
        """Impact algorithm
        """

        # Extract data
        H = get_hazard_layer(layers)    # Depth
        R = get_exposure_layer(layers)  # Building locations

        # Make the delta 10 times the size of the resolution.
        delta = abs(H.get_geotransform()[1]) * 10
        min_value, max_value = H.get_extrema()

        E = convert_line_to_points(R, delta)

        # Interpolate hazard level to building locations
        H = assign_hazard_values_to_exposure_data(H, E,
                                             attribute_name='flood_lev')

        # Extract relevant numerical data
        coordinates = E.get_geometry()
        depth = H.get_data()
        N = len(depth)

        # List attributes to carry forward to result layer
        attributes = E.get_attribute_names()

        #print attributes
        #print 'Number of population points', N

        # Calculate population impact
        road_impact = []
        num_classes = 10
        classes = range(num_classes)
        difference = (max_value - min_value) / num_classes

        for i in range(N):
            dep = float(depth[i]['flood_lev'])
            affected = classes[0]
            for level in classes:
                normalized_depth = dep - min_value
                level_value = level * difference
                if normalized_depth > level_value:
                    affected = level

            # Collect depth and calculated damage
            result_dict = {'AFFECTED': affected,
                           'DEPTH': dep}

            # Carry all original attributes forward
            for key in attributes:
                result_dict[key] = E.get_data(key, i)

            # Record result for this feature
            road_impact.append(result_dict)

        # Create report
        impact_summary = ('')

        # Create vector layer and return
        V = Vector(data=road_impact,
                   projection=E.get_projection(),
                   geometry=coordinates,
                   name='Estimated roads affected',
                   keywords={'impact_summary': impact_summary})
        return V

#.........这里部分代码省略.........
        #print 'MAX', A[245, 283], A_ref[245, 283]
        #print 'MAX: %.15f %.15f %.15f' %(A[245, 283], A_ref[245, 283])
        assert numpy.allclose(A[245, 283], A_ref[245, 283],
                              rtol=1.0e-15, atol=1.0e-15)

        #--------------
        # Exposure data
        #--------------
        # Read exposure input data for reference
        E_ref = read_layer(exposure_filename)

        # Upload to internal geonode
        exposure_layer = save_to_geonode(exposure_filename, user=self.user)
        exposure_name = '%s:%s' % (exposure_layer.workspace,
                                   exposure_layer.name)

        # Download data again
        E = download(INTERNAL_SERVER_URL, exposure_name, bbox)

        # Check exposure data against reference
        coordinates = E.get_geometry()
        coordinates_ref = E_ref.get_geometry()
        assert numpy.allclose(coordinates, coordinates_ref,
                              rtol=1.0e-12, atol=1.0e-12)

        attributes = E.get_data()
        attributes_ref = E_ref.get_data()
        for i, att in enumerate(attributes):
            att_ref = attributes_ref[i]
            for key in att:
                assert att[key] == att_ref[key]

        # Test safe's interpolation function
        I = assign_hazard_values_to_exposure_data(H, E, attribute_name='depth')
        icoordinates = I.get_geometry()

        I_ref = assign_hazard_values_to_exposure_data(H_ref, E_ref, attribute_name='depth')
        icoordinates_ref = I_ref.get_geometry()

        assert numpy.allclose(coordinates,
                              icoordinates,
                              rtol=1.0e-12, atol=1.0e-12)
        assert numpy.allclose(coordinates,
                              icoordinates_ref,
                              rtol=1.0e-12, atol=1.0e-12)

        iattributes = I.get_data()
        assert numpy.allclose(icoordinates, coordinates)

        N = len(icoordinates)
        assert N == 891

        # Set tolerance for single precision until issue #17 has been fixed
        # It appears that the single precision leads to larger interpolation
        # errors
        rtol_issue17 = 2.0e-3
        atol_issue17 = 1.0e-4

        # Verify interpolated values with test result
        for i in range(N):

            interpolated_depth_ref = I_ref.get_data()[i]['depth']
            interpolated_depth = iattributes[i]['depth']

            assert nanallclose(interpolated_depth,
                               interpolated_depth_ref,