1

I want to depict angled drill holes as a horizontal traces (plan view). I have a CSV file with X,Y,Azimuth and Horizontal length.

attribute table

Using the simple marker I can get a line pointing in the right direction, starting from the X,Y coordinates. But even when I set the size to the Horizontal Length attribute in map units, the line are not the right length. I am using UTM21N/WGS 84, so I assume it is not a projection issue.

style layer properties

  • Have you tried using the Size Assistant... tool from the same menu you tried to define the fields for the size? – Joseph Aug 25 '16 at 15:14
2

I was working on a similar problem yesterday and modified the script to meet your needs.

If you run this in the python console, the lines (along with end points) will be created and drawn as a new layer in the map window.


To run:

  • Load the csv as a point layer in QGIS

  • Edit the VARIABLES listed in the code to below to match the field names and output directories

  • Paste the code into the python console and run.


NOTE:
The line lengths will be drawn in map units, so you might want to use an EPSG for a geographic projection system.


I modified code from these sources:
Creating points based on distance and bearing from survey point using QGIS?

How to create a new empty vector layer programmatically?


import math
from PyQt4.QtCore import QVariant

# LIST OF VARIABLES:
lencolumn = 'Hor_Length'
azicolumn = 'Azimuth'
hubcode = 'HoleID'
pointcode = 'HoleID'
linecode = 'HoleID'
outputpt = r"C://temp//DELETEME//t_point.shp"
outputline = r"C://temp//DELETEME//t_line.shp"
epsg = '32621'

# CREATE POINT MEMORY LAYER
vl = QgsVectorLayer("Point?crs=EPSG:%s" % epsg, "temporary_points", "memory")
pr = vl.dataProvider()
# ADD FIELDS
pr.addAttributes([QgsField(pointcode, QVariant.String)])
vl.updateFields() # tell the vector layer to fetch changes from the provid

# CREATE LINE MEMORY LAYER
vlline = QgsVectorLayer("LineString?crs=EPSG:%s" % epsg, "temporary_lines", "memory")
prline = vlline.dataProvider()
# ADD FIELDS
prline.addAttributes([QgsField(linecode, QVariant.String)])
vlline.updateFields() # tell the vector layer to fetch changes from the provid


layer = iface.activeLayer()
iter = layer.getFeatures()
#FOR EACH INPUT FEATURE
for feature in iter:
    geom = feature.geometry()
    x = geom.asPoint()
    distance = feature[lencolumn]
    azimuth = feature[azicolumn]
    print "From Point: " + str(x)
    ####
    angle =   math.radians(azimuth)
    print "Angle: " +str(angle)     
    point = (x)
    dist_x, dist_y = \
        (distance * math.cos(angle), distance * math.sin(angle))
    #print dist_x, dist_y
    xfinal, yfinal = (point[0] + dist_x, point[1] + dist_y)
    print "output point:" 
    print xfinal, yfinal
    #insert point into empty layer
    # add a feature
    fet = QgsFeature()
    fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(xfinal, yfinal)))
    fet.setAttributes( [ str(feature[hubcode] )])
    #fet.setAttributes([str(feature.id())])
    pr.addFeatures([fet])
    # update layer's extent when new features have been added
    # because change of extent in provider is not propagated to the layer
    vl.updateExtents()
    itervl = vl.getFeatures()   
    # CREATE LINE BETWEEN POINTS WHICH SHARE THE SAME CODE(allows for multiple lines form one point)
    for v in itervl:
        print "hubcode:" + str(feature[hubcode])
        print "pointcode:" + str(v[pointcode])
        if str(v[pointcode]) == str(feature[hubcode]):  
            print "MATCH! :D"
            togeom = v.geometry()
            vx = togeom.asPoint()
            print "ToPoint: " + str(vx)
            fetline = QgsFeature()
            fetline.setGeometry(QgsGeometry.fromPolyline([QgsPoint(x),QgsPoint(vx)]))
            fetline.setAttributes( [ str(v[pointcode] )])
            prline.addFeatures([fetline])
            vlline.updateExtents()

# SAVE TO FILE / LOAD OUTPUTS TO MAP
QgsVectorFileWriter.writeAsVectorFormat(vl,outputpt,"utf-8",None,"ESRI Shapefile")
iface.addVectorLayer(outputpt, "hello_point", "ogr")
QgsVectorFileWriter.writeAsVectorFormat(vlline,outputline,"utf-8",None,"ESRI Shapefile")
iface.addVectorLayer(outputline, "hello_line", "ogr")

My test output looked like this: enter image description here

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.