Drawing Coordinates as Line in QGIS - PyQGIS

In this post I'll show some code for drawing a line from a list of coordinates, and saving it in a temporary vector layer.

You can implement the code below to get the coordinates from a text file, clipboard, or any user input.

Below are two functions that, called with the layer name and list of points as arguments, build a Line or Point temporary layer.

+1 / share if you like the post and the blog!

-----------

def createLineLay(name,lstP):
    vl = QgsVectorLayer("Linestring", name, "memory")
    #pr = vl.dataProvider()
    vl.startEditing()
    vl.addAttribute(QgsField("id", QVariant.Int))
    vl.updateFields()
    fet = QgsFeature(vl.pendingFields())
    lstPP = []
    for p in lstP:
        lstPP.append(QgsPoint(float(p[0]),float(p[1])))
    lstPP.append(QgsPoint(float(lstP[0][0]),float(lstP[0][1])))
    fet.setGeometry(QgsGeometry.fromPolyline(lstPP))
    fields = vl.pendingFields()
    fet.setFields( fields, True )
    fet['id'] = 0
    pr = vl.dataProvider()
    pr.addFeatures( [ fet ] )
    vl.commitChanges()
    QgsMapLayerRegistry.instance().addMapLayer(vl)    

def createLayPts(name, lstP):
    vl = QgsVectorLayer("Point", name, "memory")
    #pr = vl.dataProvider()
    vl.startEditing()
    vl.addAttribute(QgsField("ID", QVariant.Int))
    vl.addAttribute(QgsField("X", QVariant.Double))
    vl.addAttribute(QgsField("Y", QVariant.Double))
    vl.updateFields()
    idP = 0
    for p in lstP:
        fet = QgsFeature(vl.pendingFields())
        fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(float(p[0]),float(p[1]))))
        fields = vl.pendingFields()
        fet.setFields( fields, True )
        fet["ID"] = idP
        fet['X'] = p[0]
        fet['Y'] = p[1]
        pr = vl.dataProvider()
        pr.addFeatures( [ fet ] )
        vl.commitChanges()
        idP +=1
    QgsMapLayerRegistry.instance().addMapLayer(vl)

How to know if the polyline direction is clockwise or counter-clockwise?

I found myself puzzled with this question today, when dealing with closed lines today. I thought I could deal with this summing the horizontal deflections, but not every polygon works with this solution.
Then I found this interesting answer: 

http://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order

The answer: you have to sum over the edges, (x2 − x1)*(y2 + y1). If the result is positive the curve is clockwise, if it's negative the curve is counter-clockwise. (The result is twice the enclosed area, with a +/- convention.)

If we have a list of tuples (x,y) representing point coordinates, named lstP:

lstEdge = []
for i in range(len(lstP)-1):
    lstEdge.append((lstP[i+1][0]-lstP[i][0])*(lstP[i+1][1]+lstP[i][1]))
result = sum(lstEdge)
if result>0:
   print "Clockwise"
else:
   print "Counter-Clockwise"

    

Solving Manning's equation with Python

Manning's equation is a very common formula used in hydraulic engineering. It is an empirical formula that estimates the average velocity of open channel flow, based on a roughness coefficient.



The problem of solving Manning's formula is that it is an implicit formula, for the water depth variable, that defines the R (Hydraulic Radius).

A common approach for solving this equation is using numerical methods, as the Newton-Raphson method. An implementation for the resolution of the Manning Formula for pipes, by the Newton-Raphson approximation, is shown below.

def manningTub(depth, args):
    Q = args[0]
    diamMM = args[1]
    nMann = args[2]
    slope = args[3]
    diamM = diamMM / 1000.0
    angle = 2*acos(1-2*depth/diamM)
    area = diamM*diamM/8.0*(angle-sin(angle))
    radius = diamM/4*(1-sin(angle)/angle)
    mannR = (Q*nMann/slope**0.5)-(area*radius**(2.0/3.0))
    return mannR

def solve(f, x0, h, args):
    lastX = x0
    nextX = lastX + 10* h  # "different than lastX so loop starts OK
    while (abs(lastX - nextX) > h):  # this is how you terminate the loop - note use of abs()
        newY = f(nextX,args)  # just for debug... see what happens
        print "f(", nextX, ") = ", newY     # print out progress... again just debug
        lastX = nextX
        nextX = lastX - newY / derivative(f, lastX, h, args)  # update estimate using N-R
    return nextX

def derivative(f, x, h, args):
      return (f(x+h,args) - f(x-h,args)) / (2.0*h)  # might want to return a small non-zero if ==0

#### EXAMPLE CALL

xFound = solve(manningTub, initDepth, 0.001, [flowRate,diam,nMann,enSlope])    # call the solver

Time-Series in Python

Dealing with timeseries is a very common task in Hydrology.

One of the possibilities to process timeseries in python is to use a simple list.

For example, we can have a list of lists like this:

series1 = [ ['01/01/1900',0.0],['01/02/1900',0.1],['01/03/1900',0.3],['01/04/1900',0.4],['01/05/1900',2.2]...]

In this case, the ['01/01/1900',0.0] is composed of lists with a string representing the date, and a float number representing a value.

To properly make computations with dates, including sorting and grouping, it is necessary to interpret the string as a datetime format.

import datetime



for i in series1:

    i[0]=datetime.datetime.strptime(i[0], '%m/%d/%Y')

datetime objects accepts being sorted, making possible to sort the list based on the date, for example:

series1.sort(key=lambda x: x[0])

And we can make sums or averages based on specific months or years:

#eg. List of year 1900

lst1900 = [item for item in series1 if item[0].year==1900]



#Sum of 1900's values:

sum1900 = sum[item[1] for item in series1 if item[0].year==1900]



# avg of 1900's values

avg1900 = sum1900 / float(len(lst1900))

Extracting GPS data from pictures

It is possible to extract EXIF data from pictures, including the latitude and longitude - if it is available - with a Python native library - PIL (Python Imaging Library) .

Below is an example code of how to extract this information.

The following code opens the image file with  PIL.Image.open , then read the info with ._getexif(), naming and extracting the proper information with the functions TAGS and GPSTAGS.

The following code also includes a simple function for converting the degres minutes seconds format to decimal degrees.

from PIL import Image
import PIL.Image
from PIL.ExifTags import TAGS, GPSTAGS

def get_exif_data(image):
    """Returns a dictionary from the exif data of an PIL Image item. Also converts the GPS Tags"""
    exif_data = {}
    info = image._getexif()
    if info:
        for tag, value in info.items():
            decoded = TAGS.get(tag, tag)
            if decoded == "GPSInfo":
                gps_data = {}
                for t in value:
                    sub_decoded = GPSTAGS.get(t, t)
                    gps_data[sub_decoded] = value[t]

                exif_data[decoded] = gps_data
            else:
                exif_data[decoded] = value

    return exif_data

def _get_if_exist(data, key):
    if key in data:
        return data[key]
  
    return None
 
def _convert_to_degress(value):
    """Helper function to convert the GPS coordinates stored in the EXIF to degress in float format"""
    d0 = value[0][0]
    d1 = value[0][1]
    d = float(d0) / float(d1)

    m0 = value[1][0]
    m1 = value[1][1]
    m = float(m0) / float(m1)

    s0 = value[2][0]
    s1 = value[2][1]
    s = float(s0) / float(s1)

    return d + (m / 60.0) + (s / 3600.0)

def get_lat_lon(exif_data):
    """Returns the latitude and longitude, if available, from the provided exif_data (obtained through get_exif_data above)"""
    lat = None
    lon = None

    if "GPSInfo" in exif_data:  
        gps_info = exif_data["GPSInfo"]

        gps_latitude = _get_if_exist(gps_info, "GPSLatitude")
        gps_latitude_ref = _get_if_exist(gps_info, 'GPSLatitudeRef')
        gps_longitude = _get_if_exist(gps_info, 'GPSLongitude')
        gps_longitude_ref = _get_if_exist(gps_info, 'GPSLongitudeRef')

        if gps_latitude and gps_latitude_ref and gps_longitude and gps_longitude_ref:
            lat = _convert_to_degress(gps_latitude)
            if gps_latitude_ref != "N":                     
                lat = 0 - lat

            lon = _convert_to_degress(gps_longitude)
            if gps_longitude_ref != "E":
                lon = 0 - lon

    return lat, lon

###EXAMPLE CODE
filename = '//path//to//image.jpg'
image = PIL.Image.open(filename)
exif_data = get_exif_data(image)
latlong = get_lat_lon(exif_data)

Good and Pythonic way to iterate over files in a folder

First, we get the path by some way, most of cases with an open folder dialog.

Then we make a list with all files in that folder.

This can be useful for batch analyze / read/ write files.

import os


path = 'c://temp'
onlyfiles = [ f for f in os.listdir(path) if os.path.isfile(os.path.join(path,f))]