5

I have a point shapefile that has data for application rate, distance travelled, and swath width. In another software package it can show the points as shown.

enter image description here

I know how to make QGIS show the point data this way. However, I want to be able to use the distance travelled and swath information to create a polygon surrounding each point and assign the rate to that polygon. In the end I want a map that looks like this:

enter image description here

My other software package will display this map, but there is no way for me to export it as a shapefile so I can use it in another program. Does anyone know of a way I can get QGIS to create this type of map? I've searched around but couldn't find anything that matches what I'm trying to do.

I've worked with the software company that created the map and they do not have a way currently to export that layer. I'll see if I can get them to add that in the future, however it of course doesn't help me right now.

I've loaded a sample file with the data and information I want to use at:

Sample_Data

The data in the table I want to use is: "Swth_Wdth_" and "Distance_f" (both values are in feet) which create the rectangle around the point. I want to assign the rest of the attributes for each point to the polygon that is created (however "Rt_Apd_Liq" is the most important). If it helps, "Track_deg_" is the heading value (direction) when the point was recorded.

4
  • 2
    If you add which is your "other software package", maybe someone can also help you with how to export the data in a way QGIS can read it. Commented Apr 9, 2017 at 6:54
  • I would try to group the points by some attribute, create lines from the points and then buffer the lines with flat ends of the buffers
    – Bera
    Commented Apr 9, 2017 at 8:34
  • JohnDeere combine by any chance? My problem precisely
    – tony gil
    Commented May 26, 2018 at 17:26
  • Actually, it is from a Raven display for a manure application. It is my understanding that John Deere will export a traditional polygon shape file for its data type. Anyway, that is what I always get when I work with the data. How are you getting your data? The code below did work okay. However, it doesn't get the measurements exactly right. I think I need to convert the data to a UTM projection and try again. I ended up just interpolating the points and calling it good. Commented May 27, 2018 at 19:40

2 Answers 2

9

I generated a rectangle for each point in data according to your specification (rectangles of height = distance_f; width = swath; heading according to Track_deg). I got this result: enter image description here

There is a problem with units, so feel free to change magic constants in width and height to achieve desired result.

How to reproduce result

  1. Open your .shp layer and make sure that is the only opened layer.
  2. Open Python console. Find button to "Open Editor" (it's easier to paste code into editor than to interpreter).
  3. Paste code from below and "Run script". After a while there should be a new layer named "poly", which will contain rectangles (each has attribute Rt_Apd_Liq).

Code:

import math
from PyQt4.QtCore import QVariant

def rotate(point, angle):
    angle *= math.pi / 180.
    px, py = point
    qx = math.cos(angle) * px - math.sin(angle) * py
    qy = math.sin(angle) * px + math.cos(angle) * py
    return qx, qy

def translate(point, translation):
    return point[0] + translation[0], point[1] + translation[1]

def rectangle(width, height, angle, center):
    p1 = (-width/2, -height/2)
    p2 = (-width/2, +height/2)
    p3 = (+width/2, +height/2)
    p4 = (+width/2, -height/2)
    rotated = [rotate(p, angle) for p in [p1, p2, p3, p4]]
    translated = [translate(p, center) for p in rotated]
    points = [QgsPoint(x, y) for x, y in translated]
    return points

in_layer = QgsMapLayerRegistry.instance().mapLayers().values()[0]

layer = QgsVectorLayer('Polygon?crs=epsg:4326', 'poly' , "memory")
pr = layer.dataProvider()
pr.addAttributes([QgsField("rate", QVariant.Double)])
layer.updateFields()
for feature in in_layer.getFeatures():
    geom = feature.geometry()
    p = geom.asPoint()

    width = 0.8*0.3048 * feature.attribute('Swth_Wdth_') / 100000 # MAGIC
    length = 1.2*0.3048 * feature.attribute('Distance_f') / 100000 # MAGIC
    direction = feature.attribute('Track_deg_')
    rate = feature.attribute('Rt_Apd_Liq')
    center = (p.x(), p.y())

    poly = QgsFeature()
    points = rectangle(width, length, -direction, center)
    poly.setGeometry(QgsGeometry.fromPolygon([points]))
    poly.setAttributes([rate])
    pr.addFeatures([poly])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])
2
  • I'm not a python programmer, so I'm not sure how all of this code works. However, I did get it to work as you explained. I see the error you talk about in the width and length. Both dimensions are in feet, so I'm a little at a loss as to why the one would be (0.8*0.3048) and the other (1.2*0.3048). Shouldn't they be the same? The 0.3048 would be a conversion from feet to meters. I am assuming the 0.8 / 100000 and 1.2 /100000 has to do with converting meters to lat/long. Lat/long is not a linear conversion so depending on the direction and location those values would always be changing. Commented Apr 11, 2017 at 5:03
  • Pawel Samoraj can I kindly ask you to get an update of your excellent code so it would run in QGIS 3.10 ? I cannot get it to work... Many thanks! Commented Dec 5, 2019 at 14:31
0

If the output in your second image is a raster file, you could easily convert it to a polygon in qgis using the tool "Polygonize"

After that, to get the attributes from your point layer to the polygon layer, using "Join attributes by location"

1
  • The data in the second file is actually a bunch of polygons and not a raster. I thought about rasterizing the data, but I need to know the exact boundaries where the application happens and doesn't happen. Commented Apr 11, 2017 at 5:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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