Catchment analysis for multiple origins

mca_qgis_network_analysis

##Network_layer=vector
##Origin_layer=vector
##Network_tolerance=number 0.1
##Radius=number 2000
##Catchment_threshold=number 1000
##Output_network=output vector
##Output_catchment=output vector

# Dependencies
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

from db_manager.db_plugins.plugin import DBPlugin, Schema, Table, BaseError
from db_manager.db_plugins import createDbPlugin
from db_manager.dlg_db_error import DlgDbError
from processing.algs.qgis import postgis_utils

from shapely.ops import cascaded_union, polygonize
from shapely import geometry as sh_geom
from scipy.spatial import Delaunay
import numpy as np
import math


# Settings
#crs = QgsCoordinateReferenceSystem(27700)
otf = False
#tolerance = 0.1
#epsg = "OSGB 1936 / British National Grid"
#rad = 500

# Datasources SCAN ACTIVE LAYERS/CHOOSE CONNECTION TYPE
network_lines = processing.getObject(Network_layer)
origin_points = processing.getObject(Origin_layer)

vl = QgsVectorLayer("linestring?crs="+ crs.toWkt(), "mca_lines", "memory")
vpr = vl.dataProvider()

pl = QgsVectorLayer("polygon?crs="+ crs.toWkt(), "mca_point", "memory")
ppr = pl.dataProvider()
ppr.addAttributes([QgsField("origin", QVariant.String)]) 
pl.updateFields()


def graph_builder(network_lines, origin_points, tolerance):
 
 	# Reading crs and epsg
 	crs = network_lines.crs()
 	epsg = crs.authid()
 	
 	# Reading crs and epsg
	director = QgsLineVectorLayerDirector(network_lines, -1,'', '', '',3)
	properter = QgsDistanceArcProperter()
	director.addProperter(properter)
	builder = QgsGraphBuilder(crs,otf,tolerance,epsg)

	# Reading origins
	origins = []
	for f in origin_points.getFeatures():
		geom = f.geometry().asPoint()
		origins.append(geom)
  
	# Connect origin points to the director and build graph
	tied_origins = director.makeGraph(builder, origins)
	graph = builder.graph()


def alpha_shape(points,alpha): 
	
	# empty triangle list
	pl_lines = []
	
	# create delaunay triangulation
	pl_p_tri = Delaunay(points)

	# assess triangles
	for a, b, c in pl_p_tri.vertices:
		coord_a = pl_p_coord[a]
		coord_b = pl_p_coord[b]
		coord_c = pl_p_coord[c]
	
		# calculating length of triangle sides
		a = math.sqrt((coord_a[0] -coord_b[0]) **2 + (coord_a[1] -coord_b[1]) **2 ) 
		b = math.sqrt((coord_a[0] -coord_c[0]) **2 + (coord_a[1] -coord_c[1]) **2 ) 
		c = math.sqrt((coord_c[0] -coord_b[0]) **2 + (coord_c[1] -coord_b[1]) **2 ) 
	
		# calculating circumcircle radius
		circum_rad = (a * b * c) / math.sqrt((a+b+c) * (b+c-a) * (c+a-b) * (a+b-c))
	
		# circumcircle radius filter
		if circum_rad < alpha:
			pl_lines.append((coord_a, coord_b))
			pl_lines.append((coord_a, coord_c))
			pl_lines.append((coord_c, coord_b))
	
	# circumcircle radius filter
	pl_triangles = list(polygonize(pl_lines))
	pl_polygon = cascaded_union(pl_triangles)
	
	return pl_polygon


# Creating a graph tree per origin
i = 0

for o in tied_origins:
	
	mca_polygon_points = []
	origin_vertex_id = graph.findVertex(tied_origins[i])

    pr.addAttributes([QgsField("origin_%s" % (i+1), QVariant.Double )]) 
    vl.updateFields()
    
    (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, origin_vertex_id, 0)
    
    
	
	# Analysing the costs of the tree
    x = 0
    
    while x < len(cost):
        # lines within the radius
		if cost[x] < rad and tree[x] != -1: # lines within the radius
			outVertexId = graph.arc(tree[x]).outVertex()
			inVertexId = graph.arc(tree[x]).inVertex()
			f = QgsFeature(vl.pendingFields())
			f.setAttribute(i, cost[outVertexId]) 
			f.setGeometry(QgsGeometry.fromPolyline([graph.vertex(inVertexId).point(), graph.vertex(outVertexId).point()]))
			pr.addFeatures([f])
			
		# lines at the edge of the radius
		elif cost[x] > rad and tree[x] != -1:
			outVertexId = graph.arc(tree[x]).outVertex()
			inVertexId = graph.arc(tree[x]).inVertex()
			if cost[outVertexId] < rad:
				
				edge_line_length = rad - cost[outVertexId]
				
				edge_line_azimuth = graph.vertex(outVertexId).point().azimuth(graph.vertex(inVertexId).point())# degrees from north

				new_point_adjacent = math.sin(math.radians(edge_line_azimuth)) * edge_line_length
				new_point_opposite = math.cos(math.radians(edge_line_azimuth)) * edge_line_length

				new_point_x = graph.vertex(outVertexId).point()[0] + new_point_adjacent
				new_point_y = graph.vertex(outVertexId).point()[1] + new_point_opposite
				
				# add edge lines
				l = QgsFeature(vl.pendingFields())
				l.setAttribute(i, cost[outVertexId])
				l.setGeometry(QgsGeometry.fromPolyline([graph.vertex(outVertexId).point(),QgsPoint(new_point_x,new_point_y)]))
				pr.addFeatures([l])
				
				# add edge lines
				mca_polygon_points.append((new_point_x,new_point_y))
				
		x += 1
    
    # create concave hull
	p = QgsFeature(pl.pendingFields())
	p.setAttribute("origin", "origin_%s" % (i+1))
	p.setGeometry(QgsGeometry.fromPolygon(alpha_shape(mca_polygon_points)	
    pl.addFeatures([p])
            
    i += 1

  
vl.setCrs(crs)
vl.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(vl)
QgsMapLayerRegistry.instance().addMapLayer(pl)
print vl.isValid()
print pl.isValid()