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