Source code for safe.engine.interpolation_qgis

# coding=utf-8
"""InaSAFE Disaster risk tool by Australian Aid


.. note:: This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
     (at your option) any later version.

"""

from qgis.core import (
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsFeature,
    QgsFeatureRequest,
    QgsField,
    QgsGeometry,
    QgsRectangle,
    QgsSpatialIndex
)
from PyQt4.QtCore import QVariant

from safe.gis.qgis_vector_tools import create_layer


[docs]def interpolate_polygon_polygon(source, target, wgs84_extent): """ Transfer values from source polygon layer to the target polygon layer. This method will do a spatial join: the output layer will contain all features from target layer, with the addition of attributes of intersecting feature from the source layer. If there is not intersecting source feature, the output layer will still contain the target feature, with null attributes from the source layer. The intersection test considers only centroids of target features (not the whole polygon geometry). If more features from source layer intersect a target feature, only the first intersecting source feature will be used. :param source: Source polygon layer :type source: QgsVectorLayer :param target: Target polygon layer :type target: QgsVectorLayer :param wgs84_extent: Requested extent for analysis, in WGS84 coordinates :type wgs84_extent: QgsRectangle :return: output layer :rtype: QgsVectorLayer """ source_field_count = source.dataProvider().fields().count() target_field_count = target.dataProvider().fields().count() # Create new layer for writing resulting features. # It will contain attributes of both target and source layers result = create_layer(target) new_fields = source.dataProvider().fields().toList() new_fields.append(QgsField('polygon_id', QVariant.Int)) result.dataProvider().addAttributes(new_fields) result.updateFields() result_fields = result.dataProvider().fields() # setup transform objects between different CRS crs_wgs84 = QgsCoordinateReferenceSystem("EPSG:4326") wgs84_to_source = QgsCoordinateTransform(crs_wgs84, source.crs()) wgs84_to_target = QgsCoordinateTransform(crs_wgs84, target.crs()) source_to_target = QgsCoordinateTransform(source.crs(), target.crs()) # compute extents in CRS of layers source_extent = wgs84_to_source.transformBoundingBox(wgs84_extent) target_extent = wgs84_to_target.transformBoundingBox(wgs84_extent) # cache source layer (in CRS of target layer) source_index = QgsSpatialIndex() source_geometries = {} # key = feature ID, value = QgsGeometry source_attributes = {} for f in source.getFeatures(QgsFeatureRequest(source_extent)): f.geometry().transform(source_to_target) source_index.insertFeature(f) source_geometries[f.id()] = QgsGeometry(f.geometry()) source_attributes[f.id()] = f.attributes() # Go through all features in target layer and for each decide # whether it is intersected by any source feature result_features = [] for f in target.getFeatures(QgsFeatureRequest(target_extent)): # we use just centroids of target polygons centroid_geometry = f.geometry().centroid() centroid = centroid_geometry.asPoint() rect = QgsRectangle( centroid.x(), centroid.y(), centroid.x(), centroid.y()) ids = source_index.intersects(rect) has_matching_source = False for source_id in ids: if source_geometries[source_id].intersects(centroid_geometry): # we have found intersection between source and target f_result = QgsFeature(result_fields) f_result.setGeometry(f.geometry()) for i in xrange(target_field_count): f_result[i] = f[i] for i in xrange(source_field_count): f_result[i + target_field_count] = source_attributes[ source_id][i] f_result['polygon_id'] = source_id result_features.append(f_result) has_matching_source = True break # assuming just one source for each target feature # if there is no intersecting feature from source layer, # we will keep the source attributes null if not has_matching_source: f_result = QgsFeature(result_fields) f_result.setGeometry(f.geometry()) for i in xrange(target_field_count): f_result[i] = f[i] result_features.append(f_result) if len(result_features) == 1000: result.dataProvider().addFeatures(result_features) result_features = [] result.dataProvider().addFeatures(result_features) return result