3

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]

enter image description here

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

  • Welcome to GIS SE! As a new user please take the Tour that introduces you to this site and its protocols. – mgri Feb 14 '17 at 16:59
  • +1 Really interesting question! I would try to give you a help: please, could you upload a sample dataset (something like what I see in the image you posted)? – mgri Feb 14 '17 at 17:06
  • hey @mgri thanks! I'll add a small sample dataset to the question and the related polygon as well. – diya Feb 14 '17 at 17:29
  • Thanks! Last question: the angle of bearing is measured from the North? For example, for the selected point in the image, it should be "projected/moved" below itself? – mgri Feb 14 '17 at 17:32
  • 1
    Yep! 0/360 is North, so as you rightly pointed out, in case of the selected point it would have to projected below. I'll add that as well, thanks for the pointer! – diya Feb 14 '17 at 17:41
2

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):

enter image description here

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):

enter image description here

  • Thanks a lot @mgri!!! (even though I finally had to resort to ArcMap on a borrowed computer) But i'm really interested in the solution in QGIS.. will sift through your code while teaching myself some python! :) – diya Feb 17 '17 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.