Python QgsFeatureRequest.setFilterRect方法代码示例

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
 def writePolygonFeature(self, layer, attributes_dict):
     polygons_svg = []
     request = QgsFeatureRequest()
     request.setFilterRect(self.iface.mapCanvas().extent())
     for feature in layer.getFeatures(request):
         polygons_svg.append(self.writePolygonToSVG(feature, attributes_dict))
     return polygons_svg

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
    def __call__(self):
        if self.rect:
            rq = QgsFeatureRequest()
            rq.setFilterRect(self.rect)
            features = self.layer.getFeatures(rq)
        else:
            features = self.layer.getFeatures()

        for where in self.wheres:
            if self.DEBUG: "Has filter"
            #TODO Index lookup
#            if self.index:
#                if self.DEBUG: print "Has index"
#                min = 6163
#                max = 6164
#                iters = [iter(self.index[code]) for code in xrange(min, max + 1)]
#                features = itertools.chain(*iters)
            features = where(features)

        # TODO Clean up
        if self.limit:
            if self.DEBUG: print "Has Limit"
            for count in xrange(self.limit):
                if self.selectstatment:
                    yield self.selectstatment(features.next())
                else:
                    yield features.next()
        else:
            for f in features:
                if self.selectstatment:
                    yield self.selectstatment(f)
                else:
                    yield f

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
 def findFeaturesAt(mapPoint, layerConfig, mapTool):
     """
     To find features from a given position in a given layer
     :param mapPoint: the map position
     :param layerConfig: the layer in which we are looking for features
     :param mapTool: a QgsMapTool instance
     :return: features found in layer
     """
     if layerConfig is None:
         return None
     tolerance = layerConfig.tolerance
     if layerConfig.unit == QgsTolerance.Pixels:
         layTolerance = Finder.calcCanvasTolerance(mapTool.toCanvasCoordinates(mapPoint), layerConfig.layer, mapTool,
                                                   tolerance)
     elif layerConfig.unit == QgsTolerance.ProjectUnits:
         layTolerance = Finder.calcMapTolerance(mapPoint, layerConfig.layer, mapTool, tolerance)
     else:
         layTolerance = tolerance
     layPoint = mapTool.toLayerCoordinates(layerConfig.layer, mapPoint)
     searchRect = QgsRectangle(layPoint.x() - layTolerance, layPoint.y() - layTolerance,
                               layPoint.x() + layTolerance, layPoint.y() + layTolerance)
     request = QgsFeatureRequest()
     request.setFilterRect(searchRect)
     request.setFlags(QgsFeatureRequest.ExactIntersect)
     features = []
     for feature in layerConfig.layer.getFeatures(request):
         if layerConfig.layer.geometryType() == QGis.Polygon:
             dist, nearest, vertex = feature.geometry().closestSegmentWithContext(mapPoint)
             if QgsGeometry.fromPoint(nearest).intersects(searchRect):
                 features.append(QgsFeature(feature))
         else:
             features.append(QgsFeature(feature))
     return features

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
 def __get_text_flaeche(self, quelle, gstk, layer, fld_name):
     text = {}
     #performance! filter by bb of gstk first
     feat_req = QgsFeatureRequest()
     feat_req.setFilterRect(gstk.geometry().boundingBox())
     for feat in layer.getFeatures(feat_req):
         if feat.geometry().intersects(gstk.geometry()):
             #no fld_name defined: means yes/no only
             if fld_name is None:
                 attr_val = u'Ja'
             else:
                 attr_val = feat[fld_name]
                 #convert everything to string
                 #JSON only allows for string keys -> settingsfile
                 if isinstance( attr_val, (int, long)):
                     attr_val = unicode(attr_val)
                 elif isinstance(attr_val, float):
                     attr_val = u'{0:.0f}'.format(attr_val)
             #replace attribute values with mapping text from settings file
             if not quelle.text is None:
                 if attr_val in quelle.text:
                     attr_val = quelle.text[attr_val]
             flaeche = feat.geometry().intersection(gstk.geometry()).area()
             if fld_name in text:
                 text[attr_val] += flaeche
             else:
                 text[attr_val] = flaeche
     if len(text) < 1 and fld_name is None:
         text[u'Nein'] = 0
     elif len(text) < 1 and not fld_name is None:
         text[u'Nein'] = 0
     return text

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
def intersecting_blocks(line, destination, grid):
    """Function to fetch intersectings polygons from the grid."""
    dest_id_field = grid.fieldNameIndex('destination_id')
    request = QgsFeatureRequest()
    request.setFilterRect(line.boundingBox())
    request.setFilterExpression('"destination_id" is None')
    for feature in grid.getFeatures(request):
        if feature.geometry().intersects(line):
            grid.changeAttributeValue(feature.id(), dest_id_field, destination)

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
    def layerData(self, layer, request={}, offset=0):
        # Retrieve the data for a layer
        first = True
        data = {}
        fields = []
        fieldTypes = []
        fr = QgsFeatureRequest()
        if request:
            if 'exact' in request and request['exact']:
                fr.setFlags(QgsFeatureRequest.ExactIntersect)
            if 'nogeom' in request and request['nogeom']:
                fr.setFlags(QgsFeatureRequest.NoGeometry)
            if 'fid' in request:
                fr.setFilterFid(request['fid'])
            elif 'extents' in request:
                fr.setFilterRect(QgsRectangle(*request['extents']))
            if 'attributes' in request:
                fr.setSubsetOfAttributes(request['attributes'])

        # IMPORTANT - we do not use `for f in layer.getFeatures(fr):` as we need
        # to verify that existing attributes and geometry are correctly cleared
        # from the feature when calling nextFeature()
        it = layer.getFeatures(fr)
        f = QgsFeature()
        while it.nextFeature(f):
            if first:
                first = False
                for field in f.fields():
                    fields.append(str(field.name()))
                    fieldTypes.append(str(field.typeName()))
            if sys.version_info.major == 2:
                fielddata = dict((name, str(f[name])) for name in fields)
            else:
                fielddata = dict((name, str(f[name])) for name in fields)
            g = f.geometry()
            if not g.isEmpty():
                fielddata[geomkey] = str(g.exportToWkt())
            else:
                fielddata[geomkey] = "None"

            fielddata[fidkey] = f.id()
            id = fielddata[fields[0]]
            description = fielddata[fields[1]]
            fielddata['id'] = id
            fielddata['description'] = description
            data[f.id() + offset] = fielddata

        if 'id' not in fields:
            fields.insert(0, 'id')
        if 'description' not in fields:
            fields.insert(1, 'description')
        fields.append(fidkey)
        fields.append(geomkey)
        return fields, fieldTypes, data

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
def get_nearby_nodes(layer, node, threshold):
    """Return all nodes that has distance less than threshold from node_id.

    The list will be divided into two groups, upstream nodes and downstream
    nodes.

    :param layer: A vector point layer.
    :type layer: QGISVectorLayer

    :param node: The point/node.
    :type node: QgsFeature

    :param threshold: Distance threshold.
    :type threshold: float

    :returns: Tuple of list of nodes. (upstream_nodes, downstream_nodes).
    :rtype: tuple
    """
    id_index = layer.fieldNameIndex('id')
    node_attributes = node.attributes()
    node_id = node_attributes[id_index]
    id_index = layer.fieldNameIndex('id')
    node_type_index = layer.fieldNameIndex('node_type')
    center_node_point = node.geometry().asPoint()

    rectangle = QgsRectangle(
        center_node_point.x() - threshold,
        center_node_point.y() - threshold,
        center_node_point.x() + threshold,
        center_node_point.y() + threshold)

    # iterate through all nodes
    upstream_nodes = []
    downstream_nodes = []
    request = QgsFeatureRequest()
    request.setFilterRect(rectangle)
    for feature in layer.getFeatures(request):
        attributes = feature.attributes()
        if feature[id_index] == node_id:
            continue

        if attributes[node_type_index] == 'upstream':
            upstream_nodes.append(attributes[id_index])
        if attributes[node_type_index] == 'downstream':
            downstream_nodes.append(attributes[id_index])

    return upstream_nodes, downstream_nodes

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
def layerData(layer, request={}, offset=0):
    first = True
    data = {}
    fields = []
    fieldTypes = []
    fr = QgsFeatureRequest()
    if request:
        if 'exact' in request and request['exact']:
            fr.setFlags(QgsFeatureRequest.ExactIntersect)
        if 'nogeom' in request and request['nogeom']:
            fr.setFlags(QgsFeatureRequest.NoGeometry)
        if 'fid' in request:
            fr.setFilterFid(request['fid'])
        elif 'extents' in request:
            fr.setFilterRect(QgsRectangle(*request['extents']))
        if 'attributes' in request:
            fr.setSubsetOfAttributes(request['attributes'])

    for f in layer.getFeatures(fr):
        if first:
            first = False
            for field in f.fields():
                fields.append(str(field.name()))
                fieldTypes.append(str(field.typeName()))
        fielddata = dict((name, unicode(f[name])) for name in fields)
        g = f.geometry()
        if g:
            fielddata[geomkey] = str(g.exportToWkt())
        else:
            fielddata[geomkey] = "None"

        fielddata[fidkey] = f.id()
        id = fielddata[fields[0]]
        description = fielddata[fields[1]]
        fielddata['id'] = id
        fielddata['description'] = description
        data[f.id() + offset] = fielddata
    if 'id' not in fields:
        fields.insert(0, 'id')
    if 'description' not in fields:
        fields.insert(1, 'description')
    fields.append(fidkey)
    fields.append(geomkey)
    return fields, fieldTypes, data

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
def dissolvePolygonsOnCanvas(writer, layer):
  """dissolve polygons of the layer and clip the dissolution with base extent"""
  settings = writer.settings
  baseExtent = settings.baseExtent
  baseExtentGeom = baseExtent.geometry()
  rotation = baseExtent.rotation()
  transform = QgsCoordinateTransform(layer.crs(), settings.crs)

  combi = None
  request = QgsFeatureRequest()
  request.setFilterRect(transform.transformBoundingBox(baseExtent.boundingBox(), QgsCoordinateTransform.ReverseTransform))
  for f in layer.getFeatures(request):
    geometry = f.geometry()
    if geometry is None:
      logMessage("null geometry skipped")
      continue

    # coordinate transformation - layer crs to project crs
    geom = QgsGeometry(geometry)
    if geom.transform(transform) != 0:
      logMessage("Failed to transform geometry")
      continue

    # check if geometry intersects with the base extent (rotated rect)
    if rotation and not baseExtentGeom.intersects(geom):
      continue

    if combi:
      combi = combi.combine(geom)
    else:
      combi = geom

  # clip geom with slightly smaller extent than base extent
  # to make sure that the clipped polygon stays within the base extent
  geom = combi.intersection(baseExtent.clone().scale(0.999999).geometry())
  if geom is None:
    return None

  # check if geometry is empty
  if geom.isGeosEmpty():
    logMessage("empty geometry")
    return None

  return geom

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
    def read(self, feature_type, bbox = None, attributes = None, geometry=True, feature_filter=None):
        if not isinstance(feature_type, FeatureType):
            raise TypeError()
        lyr = self._connectlayer(feature_type)

        request = None
        if bbox or attributes is not None or not geometry or feature_filter:
            request = QgsFeatureRequest()
            if bbox:
                rect = QgsRectangle(*bbox)
                request.setFilterRect(rect)
            if attributes:
                request.setSubsetOfAttributes(attributes, lyr.pendingFields())
            if not geometry:
                request.setFlags(QgsFeatureRequest.NoGeometry)
            if feature_filter:
                request.setFilterExpression(feature_filter)
                #lyr.setSubsetString(feature_filter)

        # return listoffeatures
        # filter is maybe a QgsFeatureRequest
        # http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/vector.html#iterating-over-a-subset-of-features
        return list(lyr.getFeatures(request) if request else lyr.getFeatures())

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
    def __call__(self):
        if self.rect:
            rq = QgsFeatureRequest()
            rq.setFilterRect(self.rect)
            features = self.layer.getFeatures(rq)
        else:
            features = self.layer.getFeatures()
        
        for where in self.wheres:
            if self.DEBUG: "Has filter"
            features = where(features)

        if self.limit:
            if self.DEBUG: print "Has Limit"
            features = itertools.islice(features, 0, self.limit)
        
        if self.selectstatment:
            cols = self.selectstatment[0]
            namedcols = self.selectstatment[1]
            features = (self._project(f, *cols, **namedcols) for f in features)
        else:
            features = (self._project(f) for f in features)
        
        return features

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]

#.........这里部分代码省略.........
        raster_extent = H.dataProvider().extent()
        xmin = raster_extent.xMinimum()
        xmax = raster_extent.xMaximum()
        ymin = raster_extent.yMinimum()
        ymax = raster_extent.yMaximum()

        x_delta = (xmax - xmin) / H.width()
        x = xmin
        for i in range(H.width()):
            if abs(x - clip_xmin) < x_delta:
                # We have found the aligned raster boundary
                break
            x += x_delta
            _ = i

        y_delta = (ymax - ymin) / H.height()
        y = ymin
        for i in range(H.width()):
            if abs(y - clip_ymin) < y_delta:
                # We have found the aligned raster boundary
                break
            y += y_delta
        clip_extent = [x, y, x + width * x_delta, y + height * y_delta]

        # Clip and polygonize
        small_raster = clip_raster(
            H, width, height, QgsRectangle(*clip_extent))
        (flooded_polygon_inside, flooded_polygon_outside) = polygonize_gdal(
            small_raster, threshold_min, threshold_max)

        # Filter geometry and data using the extent
        extent = QgsRectangle(*self.extent)
        request = QgsFeatureRequest()
        request.setFilterRect(extent)

        if flooded_polygon_inside is None:
            message = tr(
                'There are no objects in the hazard layer with "value">%s.'
                'Please check the value or use other extent.' % (
                    threshold_min, ))
            raise GetDataError(message)

        #reproject the flood polygons to exposure projection
        exposure_crs = E.crs()
        exposure_authid = exposure_crs.authid()

        if hazard_authid != exposure_authid:
            flooded_polygon_inside = reproject_vector_layer(
                flooded_polygon_inside, E.crs())
            flooded_polygon_outside = reproject_vector_layer(
                flooded_polygon_outside, E.crs())

        # Clip exposure by the extent
        #extent_as_polygon = QgsGeometry().fromRect(extent)
        #no need to clip since It is using a bbox request
        #line_layer = clip_by_polygon(
        #    E,
        #    extent_as_polygon
        #)
        # Find inundated roads, mark them
        line_layer = split_by_polygon_in_out(
            E,
            flooded_polygon_inside,
            flooded_polygon_outside,
            target_field, 1, request)

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
def writeVectors(writer, legendInterface=None, progress=None):
  settings = writer.settings
  baseExtent = settings.baseExtent
  progress = progress or dummyProgress
  renderer = QgsMapRenderer()

  layers = []
  if legendInterface is None:
    for parentId in [ObjectTreeItem.ITEM_POINT, ObjectTreeItem.ITEM_LINE, ObjectTreeItem.ITEM_POLYGON]:
      for layerId, properties in settings.get(parentId, {}).iteritems():
        if properties.get("visible", False):
          layers.append([layerId, properties])
  else:
    # use vector layer order in legendInterface
    for layer in legendInterface.layers():
      if layer.type() != QgsMapLayer.VectorLayer:
        continue

      parentId = ObjectTreeItem.parentIdByLayer(layer)
      properties = settings.get(parentId, {}).get(layer.id(), {})
      if properties.get("visible", False):
        layers.append([layer.id(), properties])

  finishedLayers = 0
  for layerId, properties in layers:
    mapLayer = QgsMapLayerRegistry.instance().mapLayer(layerId)
    if mapLayer is None:
      continue

    prop = VectorPropertyReader(writer.objectTypeManager, mapLayer, properties)
    obj_mod = writer.objectTypeManager.module(prop.mod_index)
    if obj_mod is None:
      logMessage("Module not found")
      continue

    # prepare triangle mesh
    geom_type = mapLayer.geometryType()
    if geom_type == QGis.Polygon and prop.type_index == 1 and prop.isHeightRelativeToDEM():   # Overlay
      progress(None, "Initializing triangle mesh for overlay polygons")
      writer.triangleMesh()

    progress(30 + 30 * finishedLayers / len(layers), u"Writing vector layer ({0} of {1}): {2}".format(finishedLayers + 1, len(layers), mapLayer.name()))

    # write layer object
    layer = VectorLayer(writer, mapLayer, prop, obj_mod)
    writer.writeLayer(layer.layerObject(), layer.fieldNames)

    # initialize symbol rendering
    mapLayer.rendererV2().startRender(renderer.rendererContext(), mapLayer.pendingFields() if QGis.QGIS_VERSION_INT >= 20300 else mapLayer)

    # features to export
    request = QgsFeatureRequest()
    clipGeom = None
    if properties.get("radioButton_IntersectingFeatures", False):
      request.setFilterRect(layer.transform.transformBoundingBox(baseExtent.boundingBox(), QgsCoordinateTransform.ReverseTransform))
      if properties.get("checkBox_Clip"):
        extent = baseExtent.clone().scale(0.999999)   # clip with slightly smaller extent than map canvas extent
        clipGeom = extent.geometry()

    for feat in layer.features(request, clipGeom):
      # write geometry
      obj_mod.write(writer, layer, feat)   # writer.writeFeature(layer, feat, obj_mod)

      # stack attributes in writer
      if layer.writeAttrs:
        writer.addAttributes(feat.attributes())

    # write attributes
    if layer.writeAttrs:
      writer.writeAttributes()

    # write materials
    writer.writeMaterials(layer.materialManager)

    mapLayer.rendererV2().stopRender(renderer.rendererContext())
    finishedLayers += 1

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
    def run(self):
        """Experimental impact function."""
        self.validate()
        self.prepare()

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

        # Prepare Hazard Layer
        hazard_provider = self.hazard.layer.dataProvider()

        # Check affected field exists in the hazard layer
        affected_field_index = hazard_provider.fieldNameIndex(self.hazard_class_attribute)
        if affected_field_index == -1:
            message = (
                tr(
                    'Field "%s" is not present in the attribute table of the '
                    "hazard layer. Please change the Affected Field parameter in "
                    "the IF Option."
                )
                % self.hazard_class_attribute
            )
            raise GetDataError(message)

        srs = self.exposure.layer.crs().toWkt()
        exposure_provider = self.exposure.layer.dataProvider()
        exposure_fields = exposure_provider.fields()

        # Check self.exposure_class_attribute exists in exposure layer
        building_type_field_index = exposure_provider.fieldNameIndex(self.exposure_class_attribute)
        if building_type_field_index == -1:
            message = (
                tr(
                    'Field "%s" is not present in the attribute table of '
                    "the exposure layer. Please change the Building Type "
                    "Field parameter in the IF Option."
                )
                % self.exposure_class_attribute
            )
            raise GetDataError(message)

        # If target_field does not exist, add it:
        if exposure_fields.indexFromName(self.target_field) == -1:
            exposure_provider.addAttributes([QgsField(self.target_field, QVariant.Int)])
        target_field_index = exposure_provider.fieldNameIndex(self.target_field)
        exposure_fields = exposure_provider.fields()

        # Create layer to store the lines from E and extent
        building_layer = QgsVectorLayer("Polygon?crs=" + srs, "impact_buildings", "memory")
        building_provider = building_layer.dataProvider()

        # Set attributes
        building_provider.addAttributes(exposure_fields.toList())
        building_layer.startEditing()
        building_layer.commitChanges()

        # Filter geometry and data using the requested extent
        requested_extent = QgsRectangle(*self.requested_extent)

        # This is a hack - we should be setting the extent CRS
        # in the IF base class via safe/engine/core.py:calculate_impact
        # for now we assume the extent is in 4326 because it
        # is set to that from geo_extent
        # See issue #1857
        transform = QgsCoordinateTransform(
            QgsCoordinateReferenceSystem("EPSG:%i" % self._requested_extent_crs), self.hazard.layer.crs()
        )
        projected_extent = transform.transformBoundingBox(requested_extent)
        request = QgsFeatureRequest()
        request.setFilterRect(projected_extent)

        # Split building_layer by H and save as result:
        #   1) Filter from H inundated features
        #   2) Mark buildings as inundated (1) or not inundated (0)

        # make spatial index of affected polygons
        hazard_index = QgsSpatialIndex()
        hazard_geometries = {}  # key = feature id, value = geometry
        has_hazard_objects = False
        for feature in self.hazard.layer.getFeatures(request):
            value = feature[affected_field_index]
            if value not in self.hazard_class_mapping[self.wet]:
                continue
            hazard_index.insertFeature(feature)
            hazard_geometries[feature.id()] = QgsGeometry(feature.geometry())
            has_hazard_objects = True

        if not has_hazard_objects:
            message = tr(
                "There are no objects in the hazard layer with %s "
                "value in %s. Please check your data or use another "
                "attribute."
            ) % (self.hazard_class_attribute, ", ".join(self.hazard_class_mapping[self.wet]))
            raise GetDataError(message)

        features = []
        for feature in self.exposure.layer.getFeatures(request):
            building_geom = feature.geometry()
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsFeatureRequest [as 别名]
# 或者: from qgis.core.QgsFeatureRequest import setFilterRect [as 别名]
    def run(self):
        """Run the impact function.

        :returns: A new line layer with inundated roads marked.
        :type: safe_layer
        """

        target_field = self.target_field

        # Get parameters from layer's keywords
        road_class_field = self.exposure.keyword('road_class_field')
        exposure_value_mapping = self.exposure.keyword('value_mapping')

        # Get parameters from IF parameter
        threshold_min = self.parameters['min threshold'].value
        threshold_max = self.parameters['max threshold'].value

        if threshold_min > threshold_max:
            message = tr(
                'The minimal threshold is greater than the maximal specified '
                'threshold. Please check the values.')
            raise GetDataError(message)

        # reproject self.extent to the hazard projection
        hazard_crs = self.hazard.layer.crs()
        hazard_authid = hazard_crs.authid()

        if hazard_authid == 'EPSG:4326':
            viewport_extent = self.requested_extent
        else:
            geo_crs = QgsCoordinateReferenceSystem()
            geo_crs.createFromSrid(4326)
            viewport_extent = extent_to_geo_array(
                QgsRectangle(*self.requested_extent), geo_crs, hazard_crs)

        # Clip hazard raster
        small_raster = align_clip_raster(self.hazard.layer, viewport_extent)

        # Filter geometry and data using the extent
        ct = QgsCoordinateTransform(
            QgsCoordinateReferenceSystem("EPSG:4326"),
            self.exposure.layer.crs())
        extent = ct.transformBoundingBox(QgsRectangle(*self.requested_extent))
        request = QgsFeatureRequest()
        request.setFilterRect(extent)

        # create template for the output layer
        line_layer_tmp = create_layer(self.exposure.layer)
        new_field = QgsField(target_field, QVariant.Int)
        line_layer_tmp.dataProvider().addAttributes([new_field])
        line_layer_tmp.updateFields()

        # create empty output layer and load it
        filename = unique_filename(suffix='.shp')
        QgsVectorFileWriter.writeAsVectorFormat(
            line_layer_tmp, filename, "utf-8", None, "ESRI Shapefile")
        line_layer = QgsVectorLayer(filename, "flooded roads", "ogr")

        # Create vector features from the flood raster
        # For each raster cell there is one rectangular polygon
        # Data also get spatially indexed for faster operation
        index, flood_cells_map = _raster_to_vector_cells(
            small_raster,
            threshold_min,
            threshold_max,
            self.exposure.layer.crs())

        if len(flood_cells_map) == 0:
            message = tr(
                'There are no objects in the hazard layer with "value" > %s. '
                'Please check the value or use other extent.' % (
                    threshold_min, ))
            raise GetDataError(message)

        # Do the heavy work - for each road get flood polygon for that area and
        # do the intersection/difference to find out which parts are flooded
        _intersect_lines_with_vector_cells(
            self.exposure.layer,
            request,
            index,
            flood_cells_map,
            line_layer,
            target_field)

        target_field_index = line_layer.dataProvider().\
            fieldNameIndex(target_field)

        # Generate simple impact report
        epsg = get_utm_epsg(self.requested_extent[0], self.requested_extent[1])
        output_crs = QgsCoordinateReferenceSystem(epsg)
        transform = QgsCoordinateTransform(
            self.exposure.layer.crs(), output_crs)

        classes = [tr('Flooded in the threshold (m)')]
        self.init_report_var(classes)

        if line_layer.featureCount() < 1:
            raise ZeroImpactException()

        roads_data = line_layer.getFeatures()
#.........这里部分代码省略.........