Skip to content

add condition in compute_dist_to_the_hedge (windbreak_effect)

Add condition if line intersects no hedge.

def compute_dist_to_the_hedge(pnt_intersect:QgsPoint, 
                              data_hedge:QgsLineString, 
                              wind_orientation:float, 
                              hedge_id_field_name:str):
"""
Allows to calculate for each intersected point (point from `pnt_intersect`) the minimum distance to the hedge (
polyline from `data_hedge`). This function directly modifies input layer and adds field `dist`.

Parameters
---
pnt_intersect (QgsPoint): Points of the grid layer that intersect the shelter areas.  
data_hedge (QgsLineString): Layer to copy the attributes, geometry type and crs. 
wind_orientation (float): Angle in degrees indicating the origin of the wind.
hedge_id_field_name (str): Name of the id of the hedge field containing in the `data_hedge`. Defaults to "id".
Return
---
pnt_intersect (QgsPoint): Update pnt_intersect.
"""
pnt_intersect.dataProvider().addAttributes([QgsField("dist", QVariant.Double, 'double', 10, 2)])
pnt_intersect.dataProvider().addAttributes([QgsField("id", QVariant.Int)])
pnt_intersect.updateFields()
dist_idx = pnt_intersect.fields().indexFromName('dist')
id_idx = pnt_intersect.fields().indexFromName('id')
max_id = 0
attr_map = {}
wind_orientation = 90 - wind_orientation
wind_theta = math.radians(wind_orientation)
for pnt in pnt_intersect.getFeatures():
    max_id += 1
    expression = "{id} = {id_hedge}".format(id=hedge_id_field_name, id_hedge = pnt["id_hedge"])
    request = QgsFeatureRequest().setFilterExpression(expression)
    pnt_geom = pnt.geometry()
    d = 20000
    x1 = pnt_geom.asPoint().x()
    y1 = pnt_geom.asPoint().y()
    x2 = x1 + d * math.cos(wind_theta)
    y2 = y1 + d * math.sin(wind_theta)
    new_pnt = QgsPoint(x2, y2)
    new_line = QgsGeometry().fromPolyline([QgsPoint(pnt.geometry().asPoint()), new_pnt])
    hedge = next(data_hedge.getFeatures(request))
    length_dist = 0
    real_geom = []
    if new_line.intersects(hedge.geometry()):
        pnt_intersection = new_line.intersection(hedge.geometry())
        if pnt_intersection.wkbType() == QgsWkbTypes.Point:
            real_geom = pnt_intersection
        elif pnt_intersection.wkbType() == QgsWkbTypes.MultiPoint:
            real_geom = select_intersect_near_pnt(pnt_intersection, pnt_geom)
        else:
            pass
        dist = QgsGeometry().fromPolyline([QgsPoint(pnt.geometry().asPoint()), QgsPoint(real_geom.asPoint())])
        length_dist = dist.length()
        attr_map[pnt.id()] = {dist_idx: round(length_dist, 2), id_idx: max_id}
pnt_intersect.dataProvider().changeAttributeValues(attr_map)
pnt_intersect.updateFields()
return pnt_intersect