4

I have an input polygon vector layers which is essentially a single KML file containing different areas over Europe. I also have an input point vector layer which is actually flight data, which contains Lat/Lon and time stamp (24hr worth of data) of aircrafts.

What I am trying to do is to get the maximum number of aircrafts that will be inside each of these polygons over the 24hr window.

I have tried Points in Polygon, but that gives me the total number of points registered inside each polygon, but they may occur at different times of the day. What I basically want is a max of max.

EDIT

@Henhuy when i include: print(max_points, max_id, feature.id()) in the loop, i get the max values printed as follows:

enter image description here

But this is how it is shown in the attribute table:

enter image description here

In my case, row 884 (object id 885) should have a max_point of 4, but it is showing this as NULL.

  • 1
    Is there a reference in the point data which relates to an aircraft? – Knightshound Nov 18 '16 at 16:26
  • Hi, i don't fully get what you mean by the "reference", but for each time instant, i have the lat/lon of each aircraft and the aircraft registration number. However the same aircraft registration will register many other points over the course of the 24hr time window. Does that answer your question? – Amir Nov 21 '16 at 10:50
3

I tried it successfully with following code. It counts maximum points per polygon per day. For my code, point layer needs attribute "date"; then, for each date occurrences are count in dictionary and stored in layer for every polygon at end of loop. Note: New column is added to polygon layer called 'max_points' and filled with results

Code is this (I added it to python script editor in QGis, loaded both layers via selecting each layer and calling iface.activeLayer() then called the function with both layers:

Updated Code (to store values in layer):

def max_points_in_polygon_per_day(poly_layer, points_layer):
    from PyQt4.QtCore import QVariant
    poly_layer.startEditing()
    newField = QgsField('max_points', QVariant.Int)
    poly_layer.dataProvider().addAttributes([newField])
    poly_layer.updateFields()
    max_id = poly_layer.getFeatures().next().fieldNameIndex('max_points')
    for feature in poly_layer.getFeatures():
        geom = feature.geometry()
        days = {}
        for point_feature in points_layer.getFeatures():
            point_geom = point_feature.geometry()
            if point_geom is None:
                continue
            if point_geom.within(geom):
                date = point_feature.attribute('date')
                if date in days:
                    days[date] += 1
                else:
                    days[date] = 1
        if len(days) == 0:
            max_points = 0
        else:
            max_points = max(days.values())
        poly_layer.changeAttributeValue(feature.id(), max_id, max_points)
    poly_layer.commitChanges()

Looping through both feature lists could take a will and maybe is not the fastest solution - feel free to optimize my code!

EDIT

Here a short explanation how to get it to work.

First click on Python-Console:

First click on Python-Console

Then open python editor:

Then open python editor

Insert code there and click on "Run script" to reload script. Afterwards click on point layer and type following into python console:

point_layer = iface.activeLayer()

Then click on polygon layer and type:

poly_layer = iface.activeLayer()

Then type:

max_points_in_polygon_per_day(poly_layer, point_layer)

And it should print out something useful...

  • Thanks for your response Henhuy. I couldn't get this to work (being a newbie). Is there any other way to do this instead of Phyton editor? – Amir Nov 21 '16 at 13:56
  • I added a short explanation to show you how to run python script – Henhuy Nov 22 '16 at 12:46
  • Great, thanks for the detailed explanation. this is exactly what i was looking for. Just one more question: how can i add these "max" values as a new column to the attribute table of the polygon? – Amir Nov 22 '16 at 16:34
  • See new code, which stores calculated values in polygon layer. Note: Does not create new layer, but adds column to existing one! – Henhuy Nov 23 '16 at 8:49
  • Thanks. I tried the new code, which does create a new column called "max_points", but only returns "NULL" in every single row. Any thoughts? – Amir Nov 23 '16 at 10:58
2

Using the plugin 'Points2One' or 'PointsToPath' you can create lines from your points by grouping the features together on the aircraft registration number. The end product should be each flight path.

Next intersect the lines over the polygons and finally extract the centroids from the intersected line segments. Using 'Count Points in Polygon' will give you the total number of aircraft within each polygon.

  • Thanks Knightshound. I downloaded Points2One, but when i run it i get a AttributeError: 'Engine' object has no attribute 'name'. Any thoughts? – Amir Nov 21 '16 at 14:00
  • I'm afraid I can't help you there. Do you perchance have a column called 'name'? alternatively a quick search in the plugins repository shows 'PointsToPaths' which sounds like it does the same thing. Try that plugin instead – Knightshound Nov 21 '16 at 14:45
  • Thanks. I don't have a column called 'name'. i think the problem was with the format of my time-stamps which Points2One didn't seem to like. I used a dummy set and tried your approach, but using Points in Polygon, i get the number of intersection points inside each polygon. I will try playing with it a bit more. – Amir Nov 22 '16 at 15:51

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.