Python in Hydrology and Hydraulics: points

Thursday, August 31, 2017

Calculate the centroid of a polygon with python

In this post I will show a way to calculate the centroid of a non-self-intersecting closed polygon.

I used the following formulas,

$C_{\mathrm {x} }={\frac {1}{6A}}\sum _{i=0}^{n-1}(x_{i}+x_{i+1})(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})$

$C_{\mathrm {y} }={\frac {1}{6A}}\sum _{i=0}^{n-1}(y_{i}+y_{i+1})(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})$

$A={\frac {1}{2}}\sum _{i=0}^{n-1}(x_{i}\ y_{i+1}-x_{i+1}\ y_{i})\;$

as shown in https://en.wikipedia.org/wiki/Centroid#Centroid_of_a_polygon .

In the following code, the function receives the coordinates as a list of lists or as a list of tuples.

from math import sqrt

def centroid(lstP):
    sumCx = 0
    sumCy = 0
    sumAc= 0
    for i in range(len(lstP)-1):
        cX = (lstP[i][0]+lstP[i+1][0])*(lstP[i][0]*lstP[i+1][1]-lstP[i+1][0]*lstP[i][1])
        cY = (lstP[i][1]+lstP[i+1][1])*(lstP[i][0]*lstP[i+1][1]-lstP[i+1][0]*lstP[i][1])
        pA = (lstP[i][0]*lstP[i+1][1])-(lstP[i+1][0]*lstP[i][1])
        sumCx+=cX
        sumCy+=cY
        sumAc+=pA
        print cX,cY,pA
    ar = sumAc/2.0
    print ar
    centr = ((1.0/(6.0*ar))*sumCx,(1.0/(6.0*ar))*sumCy)
    return centr

Thursday, November 17, 2016

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)