本文整理汇总了Python中qgis.core.QgsFeatureRequest.setFilterRect方法的典型用法代码示例。如果您正苦于以下问题:Python QgsFeatureRequest.setFilterRect方法的具体用法?Python QgsFeatureRequest.setFilterRect怎么用?Python QgsFeatureRequest.setFilterRect使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类qgis.core.QgsFeatureRequest
的用法示例。
在下文中一共展示了QgsFeatureRequest.setFilterRect方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: writePolygonFeature
# 需要导入模块: 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
示例2: __call__
# 需要导入模块: 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
示例3: findFeaturesAt
# 需要导入模块: 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
示例4: __get_text_flaeche
# 需要导入模块: 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
示例5: intersecting_blocks
# 需要导入模块: 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)
示例6: layerData
# 需要导入模块: 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
示例7: get_nearby_nodes
# 需要导入模块: 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
示例8: layerData
# 需要导入模块: 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
示例9: dissolvePolygonsOnCanvas
# 需要导入模块: 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
示例10: read
# 需要导入模块: 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())
示例11: __call__
# 需要导入模块: 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
示例12: run
# 需要导入模块: 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)
示例13: writeVectors
# 需要导入模块: 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
示例14: run
# 需要导入模块: 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()
#.........这里部分代码省略.........
示例15: run
# 需要导入模块: 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()
#.........这里部分代码省略.........