Connecting points to polygon edges using angular bearing

Question

I am trying to map some caves located along a coastline (using QGIS 2.18.3 with GRASS 7.2.0). The survey took place on a boat, so the GPS points were taken at sea.

I have information on the angular bearing (ranging from 0-360 degrees, where 0 is North) of the caves from the boat and the lat-long position (in decimal degrees) of the boat. I plotted this information as shown below - [purple is land, yellow is position of the boat, and "rumo" is the angular bearing]

I need to find the corresponding point on the coastline based on the angular bearing and distance to shore (edge of polygon). I have tried looking at various options to do this, but as far as I can see there are ways to find the nearest distance to polygon, but not along any angular direction. The Azimuth and Distance plugin also seemed promising, but from what I understand I have to manually enter each angular bearing. Since I have >1000 such points, it is not very feasible! Seems like there is a tool in ArcMap to do this - Bearing Distance to Line. But there doesn't seem to be a direct equivalent in QGIS.

I am far from an expert in QGIS, but could be pushed into writing some code.

A sample dataset and related polygon can be found here

Answer 1

I propose a solution using PyQGIS. Since we are dealing with bearings and directions, I preferred to preliminarily work on a Projected Coordinate System and then reproject the results in a Geographic Coordinate System.

For your case, I think that this projected CRS should be correct:

Azores Central 1995 / UTM zone 26N (EPSG 3063)

So, the first task is converting your two input shapefiles to this CRS (simply right click on the layer, set the above CRS and select "Save As").

Once you have done this, you only need to run the following code as a new script from Processing:

##Points=vector point
##Polygons=vector polygon
##Search_distance=number 500

from qgis.core import *
import math

p_layer = processing.getObject(Points)
crs = p_layer.crs().toWkt()
poly_layer = processing.getObject(Polygons)

# Create the output layer
outLayer = QgsVectorLayer('Point?crs='+ crs, 'land' , 'memory')
prov = outLayer.dataProvider()
fields = p_layer.pendingFields()
prov.addAttributes(fields)
outLayer.updateFields()

all_polygons = {}
index = QgsSpatialIndex()
for ft in poly_layer.getFeatures():
    index.insertFeature(ft)
    all_polygons[ft.id()] = ft

for feat in p_layer.getFeatures():
    attrs = feat.attributes()
    point = feat.geometry().asPoint()
    points = [feat.geometry().asPoint()]
    angle = math.radians(90-(feat["bearing"]))
    dist_x, dist_y = (Search_distance * math.cos(angle), Search_distance * math.sin(angle))
    tmp_point = QgsPoint(point[0] + dist_x, point[1] + dist_y)
    points.append(tmp_point)
    tmp = QgsFeature()
    tmp.setGeometry(QgsGeometry.fromPolyline(points))
    line_geom = tmp.geometry()
    idsList = index.intersects(line_geom.boundingBox())

    for id in idsList:
        tmp_feat = all_polygons[id]
        if line_geom.intersects(tmp_feat.geometry()):
            diff = line_geom.difference(all_polygons[id].geometry())
            t_test= diff.asMultiPolyline()
            if t_test: # it's a MultiLineString
                final_point = QgsPoint(t_test[0][1][0],t_test[0][1][1])
            else: # it's a LineSring
                t_test=diff.asPolyline()
                final_point = QgsPoint(t_test[1][0],t_test[1][1])

    outFeat = QgsFeature()
    outFeat.setAttributes(attrs)
    outFeat.setGeometry(QgsGeometry.fromPoint(final_point))
    prov.addFeatures([outFeat])

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)

Two remarks:

1. In the above code, remember to adapt the name for the field which stores the bearings (so, change feat["bearing"] to feat["<name_of_the_field>"]);
2. As you can see, there is an additional parameter to insert (Search_distance): ideally, it is an infinite line from the point towards the desired direction, but I set it to 500 m which should be enough. Setting distances which lead the infinite line to intersect a polygon feature for many times may lead to the creation of the new point in the wrong location (I think you don't need to set distances bigger than mine, but you are advised).

With these premises, this will be the result (the black dashed lines are not a result, but they help to understand it):

The result is still in the projected CRS, but you can simply turn back to WGS84 if you right click on the output layer, set the desired CRS and select "Save As" (the result will be still correct):