[GRASS-SVN] r52996 - in grass-addons/grass6/raster: . r.connectivity.network

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Aug 30 11:32:23 PDT 2012


Author: sbl
Date: 2012-08-30 11:32:23 -0700 (Thu, 30 Aug 2012)
New Revision: 52996

Added:
   grass-addons/grass6/raster/r.connectivity.network/
   grass-addons/grass6/raster/r.connectivity.network/Makefile
   grass-addons/grass6/raster/r.connectivity.network/description.html
   grass-addons/grass6/raster/r.connectivity.network/r.connectivity.network
Log:
A new addon for connectivity analysis based on graph theory

Added: grass-addons/grass6/raster/r.connectivity.network/Makefile
===================================================================
--- grass-addons/grass6/raster/r.connectivity.network/Makefile	                        (rev 0)
+++ grass-addons/grass6/raster/r.connectivity.network/Makefile	2012-08-30 18:32:23 UTC (rev 52996)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.connectivity.network
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

Added: grass-addons/grass6/raster/r.connectivity.network/description.html
===================================================================
--- grass-addons/grass6/raster/r.connectivity.network/description.html	                        (rev 0)
+++ grass-addons/grass6/raster/r.connectivity.network/description.html	2012-08-30 18:32:23 UTC (rev 52996)
@@ -0,0 +1,282 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS manual: r.connectivity.network</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>r.connectivity.network</b></em>  - Compute connectivity measures for a set of habitat patches based on graph-theory
+<h2>KEYWORDS</h2>
+<dl>
+raster, vector, connectivity, network
+</dl>
+
+<h2>SYNOPSIS</h2>
+<b>r.connectivity.network</b><br>
+<b>r.connectivity.network help</b><br>
+<b>r.connectivity.network</b> [-<b>ryx</b>] <b>folder</b>=<em>string</em>  [<b>cd_cutoff</b>=<em>float</em>]   [<b>connectivity_cutoff</b>=<em>float</em>]   [<b>lnbh_cutoff</b>=<em>float</em>]   [<b>convergence_threshold</b>=<em>float</em>]   [<b>cl_thres</b>=<em>integer</em>]   [<b>euler</b>=<em>float</em>]   [<b>base</b>=<em>float</em>]   [<b>exponent</b>=<em>float</em>]   [<b>kernel_plot</b>=<em>string</em>]   [<b>overview_plot</b>=<em>string</em>]   [<b>cores</b>=<em>integer</em>]   [--<b>verbose</b>]  [--<b>quiet</b>] 
+
+<h3>Flags:</h3>
+<dl>
+<dt><b>-y</b></dt>
+<dd>Calculate edge betweenness community (default is do not).</dd>
+
+<dt><b>-x</b></dt>
+<dd>Visualise negative exponential decay kernel and exit</dd>
+
+<dt><b>--verbose</b></dt>
+<dd>Verbose module output</dd>
+<dt><b>--quiet</b></dt>
+<dd>Quiet module output</dd>
+</dl>
+
+<h3>Parameters:</h3>
+<dl>
+<dt><b>folder</b>=<em>string</em></dt>
+<dd>Folder where all (non map) output from r.connectivity.* is stored (input/output)</dd>
+
+<dt><b>cd_cutoff</b>=<em>float</em></dt>
+<dd>Maximum (cost) distance used for plotting the negative exponential decay kernel</dd>
+<dd>Default: <em>0.0</em></dd>
+
+<dt><b>connectivity_cutoff</b>=<em>float</em></dt>
+<dd>Maximum cost distance for connectivity</dd>
+<dd>Default: <em>0.0</em></dd>
+
+<dt><b>lnbh_cutoff</b>=<em>float</em></dt>
+<dd>Threshold defining a locale neighborhood (neighborhood = number of times connectivity_cutoff). Used for calculating loacle edge betweenness and neighbourhood size.</dd>
+<dd>Default: <em>0.0</em></dd>
+
+<dt><b>convergence_threshold</b>=<em>float</em></dt>
+<dd>Convergence threshold for the overview plot over the graph</dd>
+<dd>Default: <em>0.05</em></dd>
+
+<dt><b>cl_thres</b>=<em>integer</em></dt>
+<dd>Number of community levels to be traced in edge betweenness community (in addition to the number of initial clusters)</dd>
+<dd>Default: <em>0</em></dd>
+
+<dt><b>euler</b>=<em>float</em></dt>
+<dd>Euler`s number (base of the negative exponential decay kernel (<em>e ^ base * exponent</em>))</dd>
+<dd>Default: <em>2.718281828459</em></dd>
+
+<dt><b>base</b>=<em>float</em></dt>
+<dd>A factor for defining the shape of the negative exponential decay kernel (<em>e ^ base * exponent</em>)</dd>
+<dd>Default: <em>-3.0</em></dd>
+
+<dt><b>exponent</b>=<em>float</em></dt>
+<dd>Exponent of the negative exponential decay kernel (<em>e ^ base * exponent</em>)</dd>
+<dd>Default: <em>-7.5</em></dd>
+
+<dt><b>kernel_plot</b>=<em>string</em></dt>
+<dd>File name for a plot of the negative exponential decay kernel (<em>e ^ base * exponent</em>) used in analysis (requires ghostscript installed)</dd>
+
+<dt><b>overview_plot</b>=<em>string</em></dt>
+<dd>File name for a plot of an overview over network characteristics, like amount of edges, amount of network components (clusters), size of the larges network component with regards to (cost) distance, along with the connectivity threshold used in analysis (based on Bunn et al. 2000, requires ghostscript installed)</dd>
+
+<dt><b>cores</b>=<em>integer</em></dt>
+<dd>Number of cores to be used for multithreading (1 means no multithreading)</dd>
+<dd>Default: <em>1</em></dd>
+</dl>
+
+<h2>DESCRIPTION</h2>
+<dl>
+Recently, graph-theory has been characterised as an efficient and useful tool for conservation planning (e.g. Bunn et al. 2000, Calabrese & Fagan 2004, Minor & Urban 2008, Zetterberg et. al. 2010). As a part of the r.connectivity.* tool-chain, r.connectivity.network  is intended to make graph-theory more easily available to conservation planning. 
+<br><br>
+r.connectivity.network is the 2nd tool of the r.connectivity.* toolchain and performs the (core) network analysis (usig the igraph-package in R)on the network data prepared with r.connectivity.distance. This network data is analysed on graph, edge and vertex level.
+<br><br>
+Connectivity measures for the graph level are: number of vertices, number of edges, number of clusters, size of the largest cluster, average cluster size, diameter, density, modularity, number of communities, community size (in number of vertices) Network characteristics are visualised in a plot showing an  overview over number of connections, number of components and and the size of the largest network component within the network with regards to cost-distance between patches.
+<br><br>
+Connectivity measures calculated for the edges are:
+Minimum-spanning-trees, directness, edge-betweenness, local edge-betweenness, edge-betweenness-community, bridges (biconnected components), tree edges (biconnected components), and potential connectors of clusters and communities.
+<br><br>
+Connectivity measures calculated for the vertices are:
+Degree-centrality, Eigenvector-centrality, closeness centrality, vertex-betweenness, local vertex-betweenness, cluster-membership, edge betweenness community  structure, community-membership, neighbourhood size, local neighbourhood size, articulation points, and articulation.
+<br><br>
+Most measured can be calculated both on a directed and/or an undirected graph.
+<br><br>
+Analysis is based on a negative exponential decay kernel (as described e.g. in Bunn et al. (2000), which the user can modify according to the dispersal characteristics of her/his species or habitat type.
+<br><br>
+!!!This script is designed to work based on the output of r.connectivity.distance!!!
+</dl>
+
+<h2>OUTPUT NOTES</h2>
+<h3>Field naming</h3>
+<dl>
+Due to SQL-standards (especially with the dbf-driver) the length of field names is limited to ten signs. This required to abbreviate the naming of the edge and vertex measures in the vector attribute table. This was (mostly) done using the following pattern for field names (of vertex and edge measures):<br>
+<dd>weight _ measure _ (sub)graph</dd>
+</dl>
+
+<h3>(Sub)Graphs</h3>
+<dl>
+In r.connectivity.network a graph is build from the network dataset produced in r.connectivity.distance. From this (entire, directed) graph a set of (sub)graphs is extracted, specified by:
+<dd>u = for the undirected graph</dd>
+<dd>d = for the graph with only direct edges</dd>
+<dd>c = for the graph with only edges shorter than cost distance threshold</dd>
+</dl>
+
+<h3>Graph measures</h3>
+<dl>
+
+<dd><b>Number of edges</b></dd>
+<dd>The number of edges is the number of connections between pairs of vertices in the network.</dd>
+
+<dd><b>Number of vertices</b></dd>
+<dd>The number of vertices is the number of polygons (patches) from the input vector map which are analysed in the r.connectivity.*.</dd>
+
+<dd><b>Cluster</b></dd>
+<dd> A cluster is a group of vertices which are connected with each other, but not with the rest of vertices in a graph (a single isolated vertex can be a cluster as well).</dd>
+
+<dd><b>Diameter</b></dd>
+<dd>The diameter of a graph is the length of the longest geodesic.</dd>
+
+<dd><b>Density</b></dd>
+<dd>The density of a graph is the ratio of the number of edges and the number of possible edges.</dd>
+
+</dl>
+
+<h3>Edge measures</h3>
+<dl>
+<dt><b>Input parameters</b></dt>
+<dd><b>con_id</b> = ID of the directed connection between two vertices (edge)</dd>
+<dd><b>con_id_u</b> = ID of the undirected connection between two vertices (edge)</dd>
+<dd><b>from_p</b> = Patch ID (patch_id) of the start-vertex of an (directed) edge</dd>
+<dd><b>from_pop</b> = Population proxy (pop_proxy) of the start-vertex of an (directed) edge</dd>
+<dd><b>to_p</b> = Patch ID (patch_id) of the end-vertex of an (directed) edge</dd>
+<dd><b>to_pop</b> = Population proxy (pop_proxy) of the end-vertex of an (directed) edge</dd>
+
+<br><br>
+
+<dt><b>Weights</b></dt>
+
+<dd><b> Cost distance (cd)</b></dd>
+<dd>Cost distance is calculated for the directed (cd) and undirected graph (cd_u). The cd weight for directed edges is exactly like measured in r.connectivity.distance, the cd weight for undirected edges is the average cost distance of the two directed edges between a pair of vertices (mean( a--b, b-a)). The cost distance weight does not take he user defined dispersal information (the population proxy for the different patches and the dispersal kernel) into account.</dd>
+
+<dd><b>Negative exponential decay kernel flow (distk)</b></dd>
+<dd>The negative exponential decay kernel is defined by the variables <em>euler</em>, <em>base</em>, and <em>exponent</em> specified by the user. It is calculated for the directed (distk) and undirected edges(distk_u) using the following formulas:</dd>
+<dd><em>distk = euler ^ ( ( basis * ( 10 ^ exponent ) )*( cd))</em></dd>
+<dd><em>distk_u = euler ^ ( ( basis * ( 10 ^ exponent ) )*( cd_u))</em></dd>
+
+<dd><b>Maximum potential flow (mf)</b></dd>
+For the maximum potential flow (mf) weight the user defined dispersal information (the population proxy and the dispersal kernel) is used to model the flow of propagules between a pair of vertices. This weight is based on the assution that propagules disperse evenly and in full amount into all sourrounding patches. Maximum potential flow is calculated using the following formula:</dd>
+<dd>mf = disk * pop_proxy</dd>
+<dd>Maximum potential flow is calculated for the directed graph as incoming (mf_i) and outgoing (mf_o) flow as well as for the undirected graph (mf_u). Unlike the cost distance weight, mf represents the "closeness" of a pair of vertices. For some (most) algorithms it had to be inverted in order to give meaningful results (mf_inv = Inverted maximum potential flow for directed edges, mf_inv_u = Inverted maximum potential flow for undirected edges).</dd>
+
+<dd><b>Competing potential flow (cf)</b></dd>
+<dd>The competing potential flow (cf) weight is based on the work by Ranius & Roberge 2011. Like in the maximum potential flow (mf) weight, the user defined dispersal information (the population proxy and the negative exponential decay kernel) is used to model the flow of propagules between a pair of vertices. But here it is assumed that the total amount of propagules is limited to the population proxy given by the user, and that this amount is distributed according to the "attractiveness" (defined by population size and cost distance) of the surrounding patches.</dd>
+<dd>Competing potential flow is calculated for the directed (cf) and undirected (cf_u) graph. Unlike the cost distance weight, cf represents the "closeness" between a pair of vertices. For some (most) algorithms it had to be inverted in order to give meaningful results (cf_inv = Inverted competing potential flow for directed edges, cf_inv_u = Inverted competing potential flow for undirected edges).</dd>
+<br><br>
+<dd>The mf and cf weights are inverted (mf_inv, cf_inv) by multiplying them with -1 followed by linear normalising to the original maximum and minimum values.</dd>
+
+<br><br>
+
+<dt><b>Output parameters</b></dt>
+
+<dd><b>Shortest connections (isshort)</b></dd>
+<dd>In r.connectivity.network all edges are classified if they represent the shortest path from the start- to the end-vertex, which is stored as a logical value (0 = FALSE or 1 = TRUE). Shortest connections are identified for every weight (isshort_cd, isshort_mf, isshort_cf). Finally the isshort attribute is set to 1 if an edge represents the shortest path for either the cd-weight, the mf-weight or the cf-weight (or 0 otherwise).</dd>
+
+<dd><b>Biconnected components (bc)</b></dd>
+<dd>If the removal of a single vertex and its adjacent edges does not disconnect the graph (increase the number of clusters) it is part of a biconnected component.</dd>
+
+<dd><b>Bridges (is_br)</b></dd>
+<dd>A bridge is an edge who's removal increases the number of clusters in a graph. Edges are classified either as bridge (1) or not (0).</dd>
+
+<dd><b>Edge betweenness (eb)</b></dd>
+<dd>The edge betweenness (eb) value represents the number of shortest paths that go through an edge, considering the shortest paths between all possible pairs of vertices. Edge betweenness is calculated with all three weights for the entire directed graph (cd_eb_u, mf_eb_u, cf_eb_u), the undirected graph with only direct edges (cd_eb_ud, mf_eb_ud, cf_eb_ud), and the undirected graph with only direct edges shorter than cost distance threshold (cd_eb_udc, mf_eb_udc, cf_eb_udc).</dd>
+
+<dd><b>Local edge betweenness (leb)</b></dd>
+<dd>The local edge betweenness (leb) value represents the number of shortest paths that go through an edge, considering the shortest paths between all possible pairs of vertices which are shorter than the user defined local neighbourhood (lnbh_threshold). Local edge betweenness is calculated with all three weights for the entire directed graph (cd_leb_u, mf_leb_u, cf_leb_u), the undirected graph with only direct edges  (cd_leb_ud, mf_leb_ud, cf_leb_ud), and the undirected graph with only direct edges shorter than cost distance threshold (cd_leb_udc, mf_leb_udc, cf_leb_udc).</dd>
+
+<dd><b>Edge betweenness community (ebc)</b></dd>
+<dd>Some groups of vertices can be densely connected with each other but only little connected to other vertices. Such groups of relative intense connected vertices are communities. One algorithm to identify communities in a graph is edge betweenness community (ebc). Edge betweenness community is calculated in an iterative loop, where, first edge betweenness is calculated, then the edge with the largest edge betweenness value is removed, edge betweenness is recalculated and so on. The algorithm produces mainly two edge measures:</dd>
+<dd>1. the order in which edges were deleted (ebc_r)</dd>
+<dd>2. the edge betweenness value of an edge at the time of removal</dd>
+
+<dd>For large graphs calculating edge betweenness community can take long time. Therefore, it is only calculated on request (y-flag) and only for one weight (cf). It is only calculated for the entire (undirected) graph.</dd>
+
+<dd><b>Minimum spanning trees (mst)</b></dd>
+<dd>A minimum spanning tree (mst) is a subgraph consisting of the minimal possible number of edges  (= number of vertices - 1) connecting all vertices to a minimal possible number of clusters, while the sum of all edge weights is minimised. If the graph consists of more than one cluster, the result of the minimum spanning tree algorithm is a so called minimum spanning forest, consisting of the minimum spanning tree for each cluster. In the context of nature conservation, the minimum spanning tree is sometimes referred to as the "backbone" of an area network (Bunn etal. 2000). In r.connectivity.network the minimum spanning tree/forest is calculated with all weights (cd, mf, cf) and for all undirected (sub)graphs (u, ud, udc) resulting in the following (logical) edge measures cd_mst_u, cd_mst_ud, cd_mst_udc, mf_mst_u, mf_mst_ud, mf_mst_udc, cf_mst_u, cf_mst_ud, cf_mst_udc, where edges are either part of the minimum spanning tree/forest (1) or not (0).</dd>
+
+<dd><b>Potential cluster connectors (cl_pc_u)</b></dd>
+<dd>A cluster is a group of vertices which are connected with each other, but not with the rest of vertices in a graph (a single isolated vertex can be a cluster as well). Potential cluster connectors are edges which are longer than the connectivity threshold, and connect clusters in the subgraph with only edges shorter than the connectivity threshold.</dd>
+
+<dd><b>Potential community connectors (cf_ebc_pc, cf_iebc_pc)</b></dd>
+<dd>Potential community connectors are edges which connect communities identified by the edge betweenness community algorithm for a user defined community level.</dd>
+</dl>
+
+<h2>Vertex measures</h2>
+<dl>
+<dt><b>Input parameters</b></dt>
+<dd><b>Patch ID (patch_id)</b></dd>
+<dd>ID of the vertex = category value (cat) of the patch vector given by the user in r.connectivity.distance.</dd>
+
+<dd><b>Population proxy (pop_proxy)</b></dd>
+<dd>The population proxy given by the user in r.connectivity.distance.</dd>
+
+<dt><b>Output parameters</b></dt>
+<dd><b>Cluster membership (cl)</b></dd>
+<dd>A cluster is a group of vertices which are connected with each other, but not with the rest of vertices in a graph (a single isolated vertex is a cluster as well). The cluster membership value (integer) is the id of the cluster a vertex belongs to.
+
+<dd><b>Community structure and membership (cf_(i)ebc_cs, cf_(i)ebc_cl)</b></dd>
+<dd>The edge betweenness community algorithm is used to analyse the community structure of a graph. It results in a hierarchical structure which describes how a graph is split up by the edge removal process into an increasing number of communities. The hierarchical community structure is stored as the (character) vertex attribute "cf_ebc_cs". The community membership value (integer) is the id of the community (identified by edge betweeenness community) a vertex belongs to (on a user defined level of community division).</dd>
+
+<dd><b>Articulation points (art_p)</b></dd>
+<dd>Articulation points (art_p) are vertices who's removal would increase the number of clusters in a graph. Articulation points are identified for the entire directed graph (art_p_u), the undirected graph with only direct edges (art_p_ud), and the undirected graph with only direct edges shorter than cost distance threshold (art_p_udc). Articulation points (art_p) are a (logical) measure of the graph structure and do not take edge weights into account. Vertices are either articulation points (1) or not (0).</dd>
+
+<dd><b>Articulation (art)</b></dd>
+<dd>The articulation value (integer) (art) is the number of new clusters which would occur when a vertex is removed. The articulation value (art) is a measure of the graph structure and does not take edge weights into account.</dd>
+
+<dd><b>Degree centrality (deg)</b></dd>
+<dd>Degree centrality of a vertex is defined as the number of vertices connected to this vertex. This (integer) measure is a measure of the graph structure and does not take edge weights into account.. Degree centrality is calculated for all three undirected graphs (deg_u, deg_ud, deg_udc).</dd>
+
+<dd><b>Eigenvector centrality (evc)</b></dd>
+<dd>Eigenvector centrality (evc) in r.connectivity.network is a variation of the eigenvector centrality algorithm provided by the igraph library (see also: Csardi & Nepusz 2006). In r.connectivity.network eigenvector centrality is the sum of incoming potential flow of a vertex.</dd>
+
+<dd>Because eigenvector centrality takes the direction of potential flows into account it is calculated only for the directed graph and the directed graph with only edges shorter than connectivity threshold, using both cf and mf weight (cf_evc_u, cf_evc_uc, mf_evc_u, mf_evc_uc). Eigenvector centrality is stored with double precision.</dd>
+
+
+<dd><b>Closeness centrality (cl)</b></dd>
+<dd>Cloness centrality represents the number of steps which are necessary to access every other vertex from a given vertex. Closeness centrality (integer) is calculated with all three weights for the entire directed graph (cd_cl_u, mf_cl_u, cf_cl_u), the undirected graph with only direct edges  (cd_cl_ud, mf_cl_ud, cf_cl_ud), and the undirected graph with only direct edges shorter than cost distance threshold (cd_cl_udc, mf_cl_udc, cf_cl_udc).</dd>
+
+<dd><b>Vertex betweenness (vb)</b></dd>
+<dd>The vertex betweenness (vb) value represents the number of shortest paths that go through a vertex (but that do not start or end in that vertex), considering the shortest paths between all possible pairs of vertices. Vertex betweenness is calculated with all three weights for the entire directed graph (cd_vb_u, mf_vb_u, cf_vb_u), the undirected graph with only direct edges  (cd_vb_ud, mf_vb_ud, cf_vb_ud), and the undirected graph with only direct edges shorter than cost distance threshold (cd_vb_udc, mf_vb_udc, cf_vb_udc).</dd>
+
+<dd><b>Local vertex betweenness (lvb)</b></dd>
+<dd>The local vertex betweenness (lvb) value represents the number of shortest paths that go through a vertex (but that do not start or end in that vertex), considering the shortest paths between all possible pairs of vertices wich are shorter than the user defined local neighbourhood (lnbh_threshold). Local vertex betweenness is calculated with all three weights for the entire directed graph (cd_vb_u, mf_vb_u, cf_vb_u), the undirected graph with only direct edges  (cd_vb_ud, mf_vb_ud, cf_vb_ud), and the undirected graph with only direct edges shorter than cost distance threshold (cd_vb_udc, mf_vb_udc, cf_vb_udc).</dd>
+
+<dd><b>Neighbourhood size (nbh_s)</b></dd>
+<dd>The neighbourhood size (nbh_s) is the number of of other vertices which can be reached from a vertex.</dd>
+
+<dd><b>Local neighbourhood size (nbh_sl)</b></dd>
+<dd>The local neighbourhood size (nbh_sl) is the number of other vertices which can be reached from a vertex along a path which is shorter than the user defined local neighbourhood.</dd>
+
+</dl>
+
+
+<h2>REFERENCE</h2>
+<dl>
+<dt><b>Bunn, A. G., Urban, D. L. & Keitt, T. H. 2000</b>: Landscape connectivity: A conservation application of graph theory. Journal of Environmental Management (2000) 59: 265-278 <a href="http://www.sciencedirect.com/science/article/pii/S0301479700903736">http://www.sciencedirect.com/science/article/pii/S0301479700903736</a></dt>
+<dt><b>Calabrese, J. M. & Fagan, W. F. 2004</b>: A comparison-shopper's guide to connectivity metrics. Front Ecol Environ 2 (10): 529-536 <a href="http://dx.doi.org/10.1890/1540-9295(2004)002[0529:ACGTCM]2.0.CO;2">http://dx.doi.org/10.1890/1540-9295(2004)002[0529:ACGTCM]2.0.CO;2</a></dt>
+<dt><b>Minor, E. S. & Urban, D. L. 2008</b>: A Graph-Theory Framework for Evaluating Landscape Connectivity and Conservation Planning. Conservation Biology 22 (2): 297-307 <a href="http://www.uic.edu/labs/minor/minor&urban2008.pdf">http://www.uic.edu/labs/minor/minor&urban2008.pdf</a></dt>
+<dt><b>Zetterberg, A.,  Mörtberg, U. M. & Balfors, B. 2010</b>: Making graph theory operational for landscape ecological assessments, planning, and design. Landscape and Urban Planning (2010) 95: 181-191 <a href="http://www.sciencedirect.com/science/article/pii/S0169204610000204">http://www.sciencedirect.com/science/article/pii/S0169204610000204</a></dt>
+<dt><b>Ranius, T. & Roberge, J.-M. 2011</b>: Effects of intensified forestry on the landscape-scale extinction risk of dead wood dependent species. Biodiversity and Conservation 20 (13): 2867-2882<a href="http://www.springerlink.com/content/ch9qtv2665h624q4/">http://www.springerlink.com/content/ch9qtv2665h624q4/</a></dt>
+<dt></dt>
+<dt><b>Csardi G. & Nepusz T. 2006</b>: The igraph software package for complex network research, InterJournal, Complex Systems 1695. 2006. <a href="http://igraph.sf.net">http://igraph.sf.net</a></dt>
+<dt><b>Csardi, G. 2012</b>: igraph: Network analysis and visualization. <a href="http://cran.r-project.org/web/packages/igraph/index.html">http://cran.r-project.org/web/packages/igraph/index.html</a></dt>
+</dl>
+
+
+<h2>SEE ALSO</h2>
+<a href="r.connectivity.distance">r.connectivity.distance</a>, <a href="r.connectivity.corridors">r.connectivity.corridors</a><br>
+
+<h2>AUTHOR</h2>
+Stefan Blumentrath, Norwegian Institute for Nature Research (NINA)
+
+<i>Last changed: $Date$</i>
+
+</body>
+</html>

Added: grass-addons/grass6/raster/r.connectivity.network/r.connectivity.network
===================================================================
--- grass-addons/grass6/raster/r.connectivity.network/r.connectivity.network	                        (rev 0)
+++ grass-addons/grass6/raster/r.connectivity.network/r.connectivity.network	2012-08-30 18:32:23 UTC (rev 52996)
@@ -0,0 +1,1455 @@
+#!/bin/sh 
+# 
+############################################################################# 
+#
+# MODULE:       r.connectivity.network 
+# AUTHOR(S):    Stefan Blumentrath <stefan dot blumentrath at nina dot no > 
+# PURPOSE:      Compute connectivity measures for a set of habitat patches
+#               based on graph-theory (usig the igraph-package in R).
+#               
+#               Recently, graph-theory has been characterised as an efficient 
+#               and useful tool for conservation planning (e.g. Bunn et al. 2000, 
+#               Calabrese & Fagan 2004, Minor & Urban 2008, Zetterberg et. al. 2010). 
+#               As a part of the r.connectivity.* tool-chain, r.connectivity.network  
+#               is intended to make graph-theory more easily available to conservation 
+#               planning. 
+#
+#               r.connectivity.network is the 2nd tool of the r.connectivity.* toolchain 
+#               and performs the (core) network analysis (usig the igraph-package in R) 
+#               on the network data prepared with r.connectivity.distance. This network 
+#               data is analysed on graph, edge and vertex level.
+#
+#               Connectivity measures for the graph level are: number of vertices, 
+#               number of edges, number of clusters, size of the largest cluster, 
+#               average cluster size, diameter, density, modularity, number of communities, 
+#               community size (in number of vertices) Network characteristics are visualised 
+#               in a plot showing an  overview over number of connections, number of components 
+#               and and the size of the largest network component within the network with regards 
+#               to cost-distance between patches.
+#
+#               Connectivity measures calculated for the edges are:
+#               Minimum-spanning-trees, directness, edge-betweenness, local edge-betweenness, 
+#               edge-betweenness-community, bridges (biconnected components), tree edges 
+#               (biconnected components), and potential connectors of clusters and communities.
+#
+#               Connectivity measures calculated for the vertices are:
+#               Degree-centrality, Eigenvector-centrality, closeness centrality, vertex-betweenness, 
+#               local vertex-betweenness, cluster-membership, edge betweenness community  structure, 
+#               community-membership, neighbourhood size, local neighbourhood size, 
+#               articulation points, and articulation.
+#
+#               Most measured can be calculated both on a directed and/or an undirected graph.
+#
+#               Analysis is based on a negative exponential decay kernel (as described e.g. 
+#               in Bunn et al. (2000), which the user can modify according to the dispersal 
+#               characteristics of her/his species or habitat type.
+#
+#               !!!This script is designed to work based on the output of r.connectivity.distance!!!
+# 
+# COPYRIGHT:    (C) 2011 by the Norwegian Institute for Nature Research (NINA)
+# 
+#               This program is free software under the GNU General Public 
+#               License (>=v2). Read the file COPYING that comes with GRASS 
+#               for details. 
+# 
+############################################################################# 
+#
+# REQUIREMENTS:
+# R with packages igraph (version 0.6-2) and nlme (for parallel processing doMC, 
+# multicore, iterators, codetools and foreach are required as well),
+# ghostscript is required for postscript-output
+#
+#%Module 
+#% description: Compute connectivity measures for a set of habitat patches based on graph-theory
+#%End 
+#
+#%option 
+#% key: folder
+#% type: string 
+#% description: Folder where all (non map) output from r.connectivity.* is stored 
+#% required : yes
+#% guisection: Input
+#%end 
+#
+#%option 
+#% key: cd_cutoff
+#% type: double
+#% description: Maximum (cost) distance used for plotting the negative exponential decay kernel
+#% guisection: Settings
+#% required: no
+#% answer: 0.0
+#%end
+# 
+#%option 
+#% key: connectivity_cutoff
+#% type: double
+#% description: Maximum cost distance for connectivity
+#% guisection: Settings
+#% required: no
+#% answer: 0.0
+#%end
+# 
+#%option 
+#% key: lnbh_cutoff
+#% type: double
+#% description: Threshold defining a locale neighborhood (neighborhood = number of times connectivity_cutoff)
+#% guisection: Settings
+#% required: no
+#% answer: 0.0
+#%end
+# 
+#%option 
+#% key: convergence_threshold
+#% type: double
+#% description: Convergence threshold for the overview plot over the graph
+#% guisection: Settings
+#% required: no
+#% answer: 0.05
+#%end
+# 
+#%option 
+#% key: cl_thres
+#% type: integer
+#% description: Number of community levels to be traced in edge betweenness community
+#% guisection: Settings
+#% required: no
+#% answer: 0
+#%end
+#
+#%flag
+#% key: y
+#% description: Calculate edge betweenness community (default is do not).
+#% guisection: Measures
+#%end
+#
+#%option 
+#% key: euler
+#% type: double
+#% description: Euler`s number (base of the negative exponential decay kernel)
+#% guisection: Kernel
+#% required: no
+#% answer: 2.718281828459
+#%end
+# 
+#%option 
+#% key: base
+#% type: double
+#% description: A factor for defining the shape of the negative exponential decay kernel (e ^ base * exponent)
+#% guisection: Kernel
+#% required: no
+#% answer: -3.0
+#%end
+# 
+#%option 
+#% key: exponent
+#% type: double
+#% description: Exponent of the negative exponential decay kernel (e ^ base * exponent)
+#% guisection: Kernel
+#% required: no
+#% answer: -7.5
+#%end
+#
+#%flag
+#% key: x
+#% description: Visualise negative exponential decay kernel and exit
+#% guisection: Output
+#%end
+#
+#%option 
+#% key: kernel_plot
+#% type: string 
+#% description: File name for a plot of the negative exponential decay kernel (e ^ base * exponent) used in analysis (requires ghostscript installed)
+#% required : no
+#% gisprompt: new_file,file,output
+#% guisection: Output
+#%end 
+#
+#%option 
+#% key: overview_plot
+#% type: string 
+#% description: File name for a plot of an overview over network characteristics (requires ghostscript installed)
+#% required : no
+#% gisprompt: new_file,file,output
+#% guisection: Output
+#%end 
+#
+#%option 
+#% key: cores
+#% type: integer
+#% description: Number of cores to be used for computation (if <= 1 no parallelisation is applied)
+#% guisection: Parallelisation
+#% required: no
+#% answer: 1
+#%end
+# 
+
+#Input
+FOLDER=$(echo "${GIS_OPT_FOLDER}" | sed 's/\/$//g')
+
+#Input variables
+CORES="${GIS_OPT_CORES}"
+CONVERGENCE_THRESHOLD="${GIS_OPT_CONVERGENCE_THRESHOLD}"
+EULER="${GIS_OPT_EULER}"
+BASE="${GIS_OPT_BASE}"
+EXPONENT="${GIS_OPT_EXPONENT}"
+CD_CUTOFF="${GIS_OPT_CD_CUTOFF}"
+CONNECTIVITY_CUTOFF="${GIS_OPT_CONNECTIVITY_CUTOFF}"
+LNBH_CUTOFF="${GIS_OPT_LNBH_CUTOFF}"
+CL_THRES="${GIS_OPT_CL_THRES}"
+
+#Visualise negative exponential decay kernel and exit
+X_FLAG=${GIS_FLAG_X}
+
+#Calculate edge betweenness community (default is do not)
+Y_FLAG=${GIS_FLAG_Y}
+
+#
+#Check if script is started from within GRASS
+if [ -z "$GISBASE" ] ; then
+    echo "You must be in GRASS GIS to run this program." 1>&2
+    exit 1
+fi
+#
+#Pass evtl. command line arguments to gui and start it 
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    exec g.parser "$0" "$@"
+fi
+#
+#Check for log-file of the r.connectivity.* tool-chain
+if [ ! -r "${FOLDER}/r_connectivity.log" ] ; then 
+	g.message -e "r_connectivity.log file could not be found or is not readable."
+	exit 1
+fi
+
+#Read variables from log-file
+. ${FOLDER}/r_connectivity.log
+
+#Parse and assign remaining variable (Output)
+KERNEL_PLOT="${GIS_OPT_KERNEL_PLOT}"
+OVERVIEW_PLOT="${GIS_OPT_OVERVIEW_PLOT}"
+EDGES=${PREFIX}_edges
+VERTICES=${PREFIX}_vertices
+NETWORK_MEASURES="${GIS_OPT_NETWORK_MEASURES}"
+
+
+g.message -v "FOLDER is ${GIS_OPT_FOLDER}"
+g.message -v "CORES is ${GIS_OPT_CORES}"
+g.message -v "CONVERGENCE_THRESHOLD is ${GIS_OPT_CONVERGENCE_THRESHOLD}"
+g.message -v "EULER is ${GIS_OPT_EULER}"
+g.message -v "BASE is ${GIS_OPT_BASE}"
+g.message -v "EXPONENT is ${GIS_OPT_EXPONENT}"
+g.message -v "CD_CUTOFF is ${GIS_OPT_CD_CUTOFF}"
+g.message -v "CONNECTIVITY_CUTOFF is ${GIS_OPT_CONNECTIVITY_CUTOFF}"
+g.message -v "LNBH_CUTOFF is ${GIS_OPT_LNBH_CUTOFF}"
+g.message -v "KERNEL_PLOT is ${GIS_OPT_KERNEL_PLOT}"
+g.message -v "OVERVIEW_PLOT is ${GIS_OPT_OVERVIEW_PLOT}"
+g.message -v "PREFIX is ${PREFIX}"
+g.message -v "EDGES is ${PREFIX}_edges"
+g.message -v "VERTICES is ${PREFIX}_vertices"
+g.message -v "CL_THRES is ${GIS_OPT_CL_THRES}"
+
+
+#Check if R is installed
+	if [ ! -x "`which R`" ] ; then
+	    g.message -e "R is required, please install R first" 
+	    exit 1
+	fi
+
+#Visualise the negative exponential decay kernel and exit (if requested)
+if [ $X_FLAG -eq 1 ] ; then
+	
+	echo "euler <- $EULER
+	basis <- $BASE
+	exponent <- $EXPONENT
+	cd_cutoff <- $CD_CUTOFF
+	matplot((0:cd_cutoff/1000), (euler^((basis*(10^exponent))*(0:cd_cutoff))), type=\"l\", xlab=c(\"Cost distance (in 1000)\"), ylab=\"Potential flow between patches\")
+	locator(n = 1, type = \"n\")" | R --interactive --slave
+	exit 0
+	
+fi
+
+#Check if ghostscript is installed
+if [ $OVERVIEW_PLOT ] ; then
+	if [ ! -x "`which ghostscript`" ] ; then
+	    g.message -e "ghostscript is required for postscript output, please install ghostscript first" 
+	    exit 1
+	fi
+fi
+if [ $KERNEL_PLOT ]; then
+	if [ ! -x "`which ghostscript`" ] ; then
+	    g.message -e "ghostscript is required for postscript output, please install ghostscript first" 
+	    exit 1
+	fi
+fi
+
+#Check if all necessary R-packages are installed properly
+if [ $CORES -le 1 ] ; then
+missing_packages=$(echo \
+'installed_packages <- .packages(all.available = TRUE)
+missing_packages <- as.character()
+if(length(grep("igraph", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package igraph is missing")}
+if(length(grep("nlme", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package nlme is missing")}
+if(length(missing_packages)>0) {
+cat(missing_packages, sep=", ")
+}' | R --no-save --no-restore --slave)
+else
+missing_packages=$(echo \
+'installed_packages <- .packages(all.available = TRUE)
+missing_packages <- as.character()
+if(length(grep("igraph", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package igraph is missing")}
+if(length(grep("nlme", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package nlme is missing")}
+if(length(grep("doMC", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package doMC is missing")}
+if(length(grep("foreach", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package foreach is missing")}
+if(length(grep("iterators", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package iterators is missing")}
+if(length(grep("codetools", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package codetools is missing")}
+if(length(grep("multicore", installed_packages))==0) {
+missing_packages <- append(missing_packages, "Package multicore is missing")}
+if(length(missing_packages)>0) {
+cat(missing_packages, sep=", ")
+}' | R --no-save --no-restore --slave)
+fi
+if [ -n "$missing_packages" ] ; then
+	g.message -e "${missing_packages}"
+	exit 1
+fi
+
+#Check igraph version
+igraph_version=$(echo \
+'if(as.integer(as.double(sub("-", "", packageDescription("igraph", fields = c("Version"))))*100)<62) cat(packageDescription("igraph", fields = c("Version")))' | R --no-save --no-restore --slave)
+if [ -n "$igraph_version" ] ; then
+	g.message -e "The required igraph version is 0.6-2 or later, the installed version is ${igraph_version}.
+	Please download and install at least igraph version 0.6-2"
+	exit 1
+fi
+
+#Update r.connectivity.log
+cat << EOF > ${FOLDER}/r_connectivity.log
+#r.connectivity.distance
+PREFIX=$PREFIX
+PATCHES=$PATCHES
+POP_PROXY="$POP_PROXY"
+COSTS=$COSTS
+CUTOFF=$CUTOFF
+BORDER_DIST="$BORDER_DIST"
+FOLDER=$FOLDER
+EUCL_DIST=$EUCL_DIST
+#r.connectivity.network
+CORES="${GIS_OPT_CORES}"
+CONVERGENCE_THRESHOLD="${GIS_OPT_CONVERGENCE_THRESHOLD}"
+EULER="${GIS_OPT_EULER}"
+BASE="${GIS_OPT_BASE}"
+EXPONENT="${GIS_OPT_EXPONENT}"
+CD_CUTOFF="${GIS_OPT_CD_CUTOFF}"
+CONNECTIVITY_CUTOFF="${GIS_OPT_CONNECTIVITY_CUTOFF}"
+LNBH_CUTOFF="${GIS_OPT_LNBH_CUTOFF}"
+CL_THRES="${GIS_OPT_CL_THRES}"
+KERNEL_PLOT="${GIS_OPT_KERNEL_PLOT}"
+OVERVIEW_PLOT="${GIS_OPT_OVERVIEW_PLOT}"
+EDGES=${PREFIX}_edges
+VERTICES=${PREFIX}_vertices
+NETWORK_MEASURES="${FOLDER}/network_measures.csv"
+X_FLAG=${GIS_FLAG_X}
+Y_FLAG==${GIS_FLAG_Y}
+EOF
+
+
+#Start processing
+#Copy R script file to temporary folder
+cat << EOF > ${FOLDER}/network.r
+########################################################################
+###Preparation
+########################################################################
+cat("Loading variables and packages")
+
+
+library(igraph)
+library(nlme)
+library(codetools)
+library(multicore)
+library(iterators)
+library(foreach)
+library(doMC)
+
+registerDoMC()
+options(cores = $CORES)
+
+#Variables
+convergence_threshold <- $CONVERGENCE_THRESHOLD
+euler <- $EULER
+basis <- $BASE
+exponent <- $EXPONENT
+remove_indirect <- $R_FLAG
+cd_cutoff <- $CD_CUTOFF
+connectivity_cutoff <- $CONNECTIVITY_CUTOFF
+lnbh_cutoff <- $LNBH_CUTOFF
+cl_thres <- $CL_THRES
+
+#Input files
+folder <- "$FOLDER"
+
+v_input_data_1 <- paste(folder, "vertices_part_1.csv", sep="/")
+v_input_data_2 <- paste(folder, "vertices_part_2.csv", sep="/")
+
+e_input_data <- paste(folder, "edges.csv", sep="/")
+
+#Output files
+decay_kernel_ps <- "$KERNEL_PLOT"
+network_overview_ps <- "$OVERVIEW_PLOT"
+
+network_measures <- paste(folder, "network_measures.csv", sep="/")
+edge_output <- paste(folder, "edge_measures.csv", sep="/")
+edge_output_vrt <- paste(folder, "edge_measures.vrt", sep="/")
+vertex_output <- paste(folder, "vertex_measures.csv", sep="/")
+vertex_output_vrt <- paste(folder, "vertex_measures.vrt", sep="/")
+
+network_output <- paste(folder, "network_measures.csv", sep="/")
+
+if (decay_kernel_ps != "") {
+cat("Exporting a plot of the negative exponential decay kernel")
+
+########################################################################
+#Plot the negative exponential decay kernel
+ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=FALSE, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5)
+postscript(decay_kernel_ps)
+matplot((0:cd_cutoff/(10^(nchar(as.integer(cd_cutoff))-2))), (euler^((basis*(10^exponent))*(0:cd_cutoff))), type="l", xlab=paste(c("Cost distance (in ", as.character(10^(nchar(as.integer(cd_cutoff))-2)), ")"), sep="", collapse=""), ylab="Potential flow between patches")
+off <- dev.off(dev.cur())
+}
+
+########################################################################
+###Construct graph
+cat("Reading and preparing input data")
+
+#Read and join vertex data
+v_1 <- read.csv(v_input_data_1, head=TRUE, sep=';')
+v_2 <- read.csv(v_input_data_2, head=TRUE, sep=';')
+v <- merge(data.frame(patch_id=v_1[,1], centroid_x=v_1[,2], centroid_y=v_1[,3]), data.frame(patch_id=v_2[,1], pop_proxy=v_2[,2]), all.x=TRUE)
+
+#Read and group edge data
+e_df <- read.csv(e_input_data, head=TRUE, sep=';')
+#e_df <- data.frame(from_p=e_pre[,1], to_p=e_pre[,2], dist=e_pre[,3])
+#e_groups <- groupedData(dist ~ from_p | from_p/to_p, data=as.data.frame(e_df), FUN=mean)
+#e_grouped <- gsummary(e_groups, mean)
+con_id_u_pre <- unlist(mclapply(1:length(e_df\$from_p), function(x) paste(sort(c(e_df\$from_p[x], e_df\$to_p[x]))[1], sort(c(e_df\$from_p[x], e_df\$to_p[x]))[2], sep="_")))
+con_id_u_pre_df <- data.frame(con_id_u_pre=unique(con_id_u_pre), con_id_u=1:length(unique(con_id_u_pre)))
+e_grouped_df_pre <- merge(data.frame(con_id_u_pre=con_id_u_pre, con_id=1:length(e_df\$from_p), from_p=e_df\$from_p, to_p=e_df\$to_p, dist=e_df\$cost.dist), con_id_u_pre_df, all.x=TRUE)
+e_groups_ud <- groupedData(dist ~ 1 | con_id_u, data=e_grouped_df_pre, FUN=mean)
+e_grouped_ud <- gsummary(e_groups_ud, mean)
+e_grouped_ud_df <- data.frame(con_id_u=as.vector(e_grouped_ud\$con_id_u), dist_ud=as.vector(e_grouped_ud\$dist))
+e_grouped_df <- merge(data.frame(con_id=as.integer(e_grouped_df_pre\$con_id), con_id_u=as.integer(e_grouped_df_pre\$con_id_u), from_p=as.integer(e_grouped_df_pre\$from_p), to_p=as.integer(e_grouped_df_pre\$to_p), dist=as.double(e_grouped_df_pre\$dist)), e_grouped_ud_df, all.x=TRUE)
+
+##Remove connections longer than cost distance threshold if requested
+#if(remove_longer_cutoff == 1) {
+#con_id_true <-  e_grouped_df\$con_id[grep(TRUE, e_grouped_df\$dist_ud<cd_cutoff)]
+#e_grouped_df_pre <- merge(data.frame(con_id=con_id_true), e_grouped_df)
+#con_id_u_df <- data.frame(con_id_u_old=unique(e_grouped_df_pre\$con_id_u), con_id_u=1:length(unique(e_grouped_df_pre\$con_id_u)))
+#e_grouped_df <- merge(data.frame(con_id_old=e_grouped_df_pre\$con_id, con_id_u_old=e_grouped_df_pre\$con_id_u, con_id=1:length(e_grouped_df_pre\$con_id), from_p=e_grouped_df_pre\$from_p, to_p=e_grouped_df_pre\$to_p, dist=e_grouped_df_pre\$dist, dist_ud=e_grouped_df_pre\$dist_ud), con_id_u_df, all.x=TRUE)
+#}
+
+
+#Merge vertex and grouped edge data
+from_pop <- merge(data.frame(con_id=e_grouped_df\$con_id, from_p=e_grouped_df\$from_p), data.frame(from_p=as.integer(v\$patch_id), from_pop=as.double(v\$pop_proxy), from_centroid_x=as.numeric(v\$centroid_x), from_centroid_y=as.numeric(v\$centroid_y)), all.x=TRUE, sort=FALSE)
+to_pop <- merge(data.frame(con_id=e_grouped_df\$con_id, to_p=e_grouped_df\$to_p), data.frame(to_p=as.integer(v\$patch_id), to_pop=as.double(v\$pop_proxy), to_centroid_x=as.numeric(v\$centroid_x), to_centroid_y=as.numeric(v\$centroid_y)), all.x=TRUE, sort=FALSE)
+e <- merge(merge(from_pop, to_pop, sort=TRUE), data.frame(con_id=e_grouped_df\$con_id, con_id_u=e_grouped_df\$con_id_u, cost_distance=e_grouped_df\$dist, cd_u=e_grouped_df\$dist_ud), by="con_id", sort=TRUE)
+
+#Calculate attributes representing proxies for potential flow between patches
+distance_weight_e <- euler^((basis*(10^exponent))*(e\$cost_distance))
+distance_weight_e_ud <- euler^((basis*(10^exponent))*(e\$cd_u))
+#out_areal_distance_weight_e <- e\$from_pop*distance_weight_e
+#in_areal_distance_weight_e <- e\$to_pop*distance_weight_e
+
+mf_o <- e\$from_pop*distance_weight_e
+mf_i <- e\$to_pop*distance_weight_e
+mf_u <- e\$from_pop*distance_weight_e_ud+e\$from_pop*distance_weight_e_ud
+
+mf_o_inv <- mf_o^-1
+mf_i_inv <- mf_i^-1
+mf_inv_u <- ((mf_u*-1)-min(mf_u*-1))*(ifelse(max(mf_u)>10000000000000000,10000000000000000,max(mf_u))-ifelse(min(mf_u)<0.000000000001,0.000000000001,min(mf_u)))/(max(mf_u*-1)-min(mf_u*-1))+ifelse(min(mf_u)<0.000000000001,0.000000000001,min(mf_u))
+
+
+flow_df <- data.frame(from_p=e\$from_p, mf_i=mf_i)
+sum_mf_i_groups <- groupedData(mf_i ~ from_p | from_p, data=as.data.frame(flow_df), FUN=sum)
+sum_mf_i_pre <- gsummary(sum_mf_i_groups, sum)
+sum_mf_i <- data.frame(from_p=as.integer(as.character(sum_mf_i_pre[,1])), sum_mf_i=as.double(sum_mf_i_pre[,2]))
+
+tmp <- merge(merge(e, data.frame(con_id=e\$con_id, distance_weight_e_ud=distance_weight_e_ud, distance_weight_e=distance_weight_e, mf_o=mf_o, mf_o_inv=mf_o_inv, mf_i=mf_i, mf_i_inv=mf_i_inv, mf_u=mf_u, mf_inv_u=mf_inv_u)), sum_mf_i)
+cf <- tmp\$mf_o*(tmp\$mf_i/tmp\$sum_mf_i)
+#cf <- (cf-rep(min(cf),length(cf)))*(1/(rep(max(cf),length(cf))-rep(min(cf),length(cf))))
+#cf_inv <- cf^-1
+#igraphs limits for weights in edge.betweenness function are aproximately max 12600000000000000 and min 0.000000000001
+cf_inv <- ((cf*-1)-min(cf*-1))*(ifelse(max(cf)>10000000000000000,10000000000000000,max(cf))-ifelse(min(cf)<0.000000000001,0.000000000001,min(cf)))/(max(cf*-1)-min(cf*-1))+ifelse(min(cf)<0.000000000001,0.000000000001,min(cf))
+
+e <- merge(tmp, data.frame(con_id=e\$con_id, cf=cf, cf_inv=cf_inv), by="con_id", sort=TRUE)
+
+cf_u_df <- data.frame(con_id_u=e\$con_id_u, cf=e\$cf, cf_inv=e\$cf_inv)
+cf_u_groups <- groupedData(cf ~ 1 | con_id_u, data=cf_u_df, FUN=sum)
+cf_u_grouped <- gsummary(cf_u_groups, sum)
+
+e <- merge(e, data.frame(con_id_u=as.vector(cf_u_grouped\$con_id_u), cf_u=as.vector(cf_u_grouped\$cf), cf_inv_u=((as.vector(cf_u_grouped\$cf)*-1)-min(as.vector(cf_u_grouped\$cf)*-1))*(ifelse(max(as.vector(cf_u_grouped\$cf))>10000000000000000,10000000000000000,max(as.vector(cf_u_grouped\$cf)))-ifelse(min(as.vector(cf_u_grouped\$cf))<0.000000000001,0.000000000001,min(as.vector(cf_u_grouped\$cf))))/(max(as.vector(cf_u_grouped\$cf)*-1)-min(as.vector(cf_u_grouped\$cf)*-1))+ifelse(min(as.vector(cf_u_grouped\$cf))<0.000000000001,0.000000000001,min(as.vector(cf_u_grouped\$cf)))), all.x=TRUE)
+
+#Collapse edge data for undirected graph
+con_id_u_unique <- unique(e\$con_id_u)
+ud_groups <- groupedData(con_id ~ 1 | con_id_u, data=data.frame(con_id_u=e\$con_id_u, con_id=e\$con_id), FUN=mean)
+ud_grouped <- gsummary(ud_groups)
+ud_grouped_df <- data.frame(con_id_u=as.vector(ud_grouped\$con_id_u), con_id_avg=as.vector(ud_grouped\$con_id))
+e_ud_pre <- merge(e, ud_grouped_df, all.x=TRUE)
+first_con_id <- data.frame(con_id=e_ud_pre\$con_id[e_ud_pre\$con_id_avg>e_ud_pre\$con_id])
+e_ud <- merge(first_con_id, e)
+
+#Merge vertices and edges to graph-object
+#Build directed graph
+cat("Building the graph")
+
+g <- graph.empty()
+g <- add.vertices(g, nrow(v), patch_id=as.character(v\$patch_id), pop_proxy=as.numeric(v\$pop_proxy), centroid_x=as.numeric(v\$centroid_x), centroid_y=as.numeric(v\$centroid_y))
+names <- V(g)\$patch_id
+ids <- 1:length(names)
+names(ids) <- names
+from <- as.character(e\$from_p)
+to <- as.character(e\$to_p)
+edges <- matrix(c(ids[from], ids[to]), nc=2)
+g <- add.edges(g, t(edges), con_id=e\$con_id, con_id_u=e\$con_id_u, from_p=e\$from_p, from_pop=e\$from_pop, from_p_centroid_x=e\$from_centroid_x, from_p_centroid_y=e\$from_centroid_y,  to_p=e\$to_p, to_pop=e\$to_pop, to_p_centroid_x=e\$to_centroid_x, to_p_centroid_y=e\$to_centroid_y, cost_distance=e\$cost_distance, cd_u=e\$cd_u, distance_weight_e=e\$distance_weight_e, distance_weight_e_ud=e\$distance_weight_e_ud, mf_o=e\$mf_o, mf_o_inv=e\$mf_o_inv, mf_i=e\$mf_i, mf_i_inv=e\$mf_i_inv, mf_u=e\$mf_u, mf_inv_u=e\$mf_inv_u, cf=e\$cf, cf_inv=e\$cf_inv, cf_u=e\$cf_u, cf_inv_u=e\$cf_inv_u)
+
+#Build undirected graph
+g_ud <- as.undirected(graph.empty())
+g_ud <- add.vertices(g_ud, nrow(v), patch_id=as.character(v\$patch_id), pop_proxy=as.numeric(v\$pop_proxy), centroid_x=as.numeric(v\$centroid_x), centroid_y=as.numeric(v\$centroid_y))
+names <- V(g_ud)\$patch_id
+ids <- 1:length(names)
+names(ids) <- names
+from <- as.character(e_ud\$from_p)
+to <- as.character(e_ud\$to_p)
+edges <- matrix(c(ids[from], ids[to]), nc=2)
+g_ud <- add.edges(g_ud, t(edges), con_id=e_ud\$con_id, con_id_u=e_ud\$con_id_u, from_p=e_ud\$from_p, from_pop=e_ud\$from_pop, from_p_centroid_x=e_ud\$from_centroid_x, from_p_centroid_y=e_ud\$from_centroid_y,  to_p=e_ud\$to_p, to_pop=e_ud\$to_pop, to_p_centroid_x=e_ud\$to_centroid_x, to_p_centroid_y=e_ud\$to_centroid_y, cost_distance=e_ud\$cost_distance, cd_u=e_ud\$cd_u, distance_weight_e=e_ud\$distance_weight_e, mf_o=e_ud\$mf_o, mf_o_inv=e_ud\$mf_o_inv, mf_i=e_ud\$mf_i, mf_i_inv=e_ud\$mf_i_inv, mf_u=e_ud\$mf_u, mf_inv_u=e_ud\$mf_inv_u, cf=e_ud\$cf, cf_inv=e_ud\$cf_inv, cf_u=e_ud\$cf_u, cf_inv_u=e_ud\$cf_inv_u)
+
+#Classify connection with regards to shortest paths
+#con_id_u <- unlist(lapply(0:(length(E(g))-1), function(x) E(g, path=c(sort(array(c(get.edge(g, x)[1], get.edge(g, x)[2])))[1], sort(array(c(get.edge(g, x)[1], get.edge(g, x)[2])))[2]))))
+#l1 <- unlist(lapply(1:length(V(g_ud)), function(x) rep(x, length(V(g_ud)))))
+#l2 <- rep(1:length(V(g_ud)), length(V(g_ud)))
+#fa <- function(a) get.shortest.paths(g_ud, l1[a], l2[a], weights=E(g_ud)\$cd_u)
+#sp <- fa(1:length(l1))
+#Set edge attributes
+#E(g)\$con_id_u <- con_id_u
+
+#E(g)\$cd_u <- unlist(mclapply(0:(length(E(g))-1), function(x) (E(g)[x]\$cost_distance+E(g, path=c(get.edge(g, x)[2], get.edge(g, x)[1]))\$cost_distance)/2))
+
+####################################
+###Can the following be vectorised????
+
+E(g_ud)\$isshort_cd <- as.integer(unlist(mclapply(1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=get.edge(g_ud, x)[1], to=get.edge(g_ud, x)[2], weights=E(g_ud)\$cd_u)[[1]])))==2)
+E(g_ud)\$isshort_mf <- as.integer(unlist(mclapply(1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=get.edge(g_ud, x)[1], to=get.edge(g_ud, x)[2], weights=E(g_ud)\$mf_inv_u)[[1]])))==2)
+E(g_ud)\$isshort_cf <- as.integer(unlist(mclapply(1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=get.edge(g_ud, x)[1], to=get.edge(g_ud, x)[2], weights=E(g_ud)\$cf_inv_u)[[1]])))==2)
+E(g_ud)\$isshort <- ifelse((E(g_ud)\$isshort_cd==0 & E(g_ud)\$isshort_mf==0 & E(g_ud)\$isshort_cf==0), 0, 1)
+####################################
+
+#Add edge attributes to directed graph (for export)
+#E(g)\$isshort[grep(TRUE, E(g)\$con_id_u %in% E(g_ud)\$con_id_u[grep(TRUE, E(g_ud)\$isshort)])] <- TRUE
+#E(g)\$isshort[grep(TRUE, E(g)\$con_id_u %in% E(g_ud)\$con_id_u[grep(FALSE, E(g_ud)\$isshort)])] <- FALSE
+
+#Remove indirect connections if requested
+#if(remove_indirect == 1) {
+g_ud_d <- delete.edges(g_ud, E(g_ud)[E(g_ud)\$isshort==FALSE])
+#g_ud <- delete.edges(g_ud, E(g_ud)[E(g_ud)\$isshort==FALSE])
+#}
+
+#Remove connections above connectivity threshold if requested
+if(connectivity_cutoff >= 0.0) {
+g_ud_cd <- delete.edges(g_ud, E(g_ud)[grep(TRUE, E(g_ud)\$cd_u>=connectivity_cutoff)])
+g_ud_d_cd <- delete.edges(g_ud_d, E(g_ud_d)[grep(TRUE, E(g_ud_d)\$cd_u>=connectivity_cutoff)])
+}
+
+########################################################################
+###Analysis on graph level
+########################################################################
+
+cat("Starting analysis on graph level")
+
+###graph measures
+vertices_n <- length(V(g_ud_d))
+edges_n <- length(E(g_ud))
+edges_n_d <- length(E(g_ud_d))
+edges_n_cd <- length(E(g_ud_cd))
+edges_n_d_cd <- length(E(g_ud_d_cd))
+
+
+########################################################################
+######Calculate number of clusters
+#On the undirected graph with only direct edges and on the undirected graph with only direct edges shorter cost distance threshold
+cl_no_d <- clusters(g_ud_d)\$no
+cl_no_d_cd <- clusters(g_ud_d_cd)\$no
+
+cls_size_d <- as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d)\$membership, pop_size=V(g_ud_d)\$pop_proxy), order.groups=TRUE, FUN=sum), sum)\$pop_size)
+cls_size_d_cd <- as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d_cd)\$membership, pop_size=V(g_ud_d_cd)\$pop_proxy), order.groups=TRUE, FUN=sum), sum)\$pop_size)
+
+cl_max_size_d <- max(cls_size_d)
+cl_max_size_d_cd <- max(cls_size_d_cd)
+
+cl_mean_size_d <- mean(cls_size_d)
+cl_mean_size_d_cd <- mean(cls_size_d_cd)
+
+diam <- diameter(g_ud, directed=FALSE, unconnected=TRUE, weights=E(g_ud)\$cd_u)
+diam_d <- diameter(g_ud_d, directed=FALSE, unconnected=TRUE, weights=E(g_ud_d)\$cd_u)
+diam_d_cd <- diameter(g_ud_d_cd, directed=FALSE, unconnected=TRUE, weights=E(g_ud_d_cd)\$cd_u)
+
+density <- graph.density(g, loops=FALSE)
+density_u <- graph.density(g_ud, loops=FALSE)
+density_ud <- graph.density(g_ud_d, loops=FALSE)
+density_udc <- graph.density(g_ud_d_cd, loops=FALSE)
+
+if (network_overview_ps != "") {
+######Edge removal operations
+idx <- sort(E(g_ud)\$cd_u, decreasing=TRUE, na.last=NA, index.return=TRUE)\$ix
+cl_del_count <- unlist(mclapply(1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no/vertices_n)>=convergence_threshold||E(g_ud)\$cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no,NA)))
+cl_del_max_size <- unlist(mclapply(1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no/vertices_n)>=convergence_threshold||E(g_ud)\$cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),max(as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$membership, pop_size=V(g_ud)\$pop_proxy), order.groups=TRUE, FUN=sum), sum)\$pop_size)),NA)))
+cl_del_diam <- unlist(mclapply(1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no/vertices_n)>=convergence_threshold||E(g_ud)\$cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),diameter(delete.edges(g_ud, E(g_ud)[idx[1:x]]), directed=FALSE, unconnected=TRUE, weights=E(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$cd_u),NA)))
+extract <- sort(grep(FALSE, is.na(cl_del_count)), decreasing=TRUE)
+df_edgeremoval <- data.frame(distance=E(g_ud)\$cd_u[((idx[extract]))], cl_del_count=cl_del_count[extract], cl_del_max_size=cl_del_max_size[extract], cl_del_diam=cl_del_diam[extract])
+dist_inv_ix <- sort(df_edgeremoval\$distance, decreasing=FALSE, na.last=NA, index.return=TRUE)\$ix
+
+ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=FALSE, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5)
+postscript(network_overview_ps)
+
+###Check axis lables!!!
+matplot(df_edgeremoval\$distance/(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), df_edgeremoval\$cl_del_count*100/vertices_n, type="l", ylab="", xlab=c("Connectivity threshold", paste("(Cost distance between patches in ", as.character(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), ")", sep="")), lty=1, yaxt="n", yaxs="r", ylim=c(0, 100), yaxs="i", xaxs="i")
+if (connectivity_cutoff>0) abline(v=connectivity_cutoff/10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2), col="red", lty=3)
+ifelse(connectivity_cutoff>0, legend("topleft", c("Clusters (in % of maximum possible clusters)", "Size of the largest cluster (in % of total population size)", "Number of edges (in % of maximum possible number of edges)", "Diameter (in % of diameter of the entire graph)", "Connectivity threshold used in analysis"), lty=c(1, 2, 3, 4, 3), col=c("black", "black", "black", "black", "red"), inset=0.005, bty="o", box.lty=0, bg="White"), legend("topleft", c("Clusters (in % of maximum possible clusters)", "Size of the largest cluster (in % of total population size)", "Number of edges (in % of maximum possible number of edges)", "Diameter (in % of diameter of the entire graph)"), lty=c(1, 2, 3, 4), col=c("black", "black", "black", "black"), inset=0.005, bty="o", box.lty=0, bg="White"))
+axis(2, seq.int(0, 100, 25), labels=c("0 %", "25 %", "50 %", "75 %", "100 %"), yaxs="r")
+matplot(df_edgeremoval\$distance/(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), df_edgeremoval\$cl_del_max_size*100/sum(V(g)\$pop_proxy), type="l", lty=2, yaxt="n", yaxs="i", add=TRUE)
+lines(df_edgeremoval\$distance/(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), dist_inv_ix*100/edges_n, type="l", lty=3)
+lines(df_edgeremoval\$distance/(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), df_edgeremoval\$cl_del_diam*100/diam_d, type="l", lty=4)
+#Lable axis 4!!!
+#axis(4, seq.int(0, 100, 25), yaxs="i", labels=seq.int(0, ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2)), ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2))/4))
+off <- dev.off(dev.cur())
+
+#idx <- sort(E(g_ud)\$cd_u, decreasing=TRUE, na.last=NA, index.return=TRUE)\$ix
+#cl_del_count <- unlist(mclapply(1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no/vertices_n)>=convergence_threshold,clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no,NA)))
+#cl_del_max_size <- unlist(mclapply(1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no/vertices_n)>=convergence_threshold,max(as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$membership, pop_size=V(g_ud)\$pop_proxy), order.groups=TRUE, FUN=sum), sum)\$pop_size)),NA)))
+#cl_del_diam <- unlist(mclapply(1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$no/vertices_n)>=convergence_threshold,diameter(delete.edges(g_ud, E(g_ud)[idx[1:x]]), directed=FALSE, unconnected=TRUE, weights=E(delete.edges(g_ud, E(g_ud)[idx[1:x]]))\$cd_u),NA)))
+#extract <- sort(grep(FALSE, is.na(cl_del_count)), decreasing=TRUE)
+#df_edgeremoval <- data.frame(distance=E(g_ud)\$cd_u[((idx[extract]))], cl_del_count=cl_del_count[extract], cl_del_max_size=cl_del_max_size[extract], cl_del_diam=cl_del_diam[extract])
+#dist_inv_ix <- sort(df_edgeremoval\$distance, decreasing=FALSE, na.last=NA, index.return=TRUE)\$ix
+
+#ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=FALSE, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5)
+#postscript(network_overview_ps)
+
+####Check axis lables!!!
+#matplot(df_edgeremoval\$distance/(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), df_edgeremoval\$cl_del_count*100/vertices_n, type="l", xlab=c("Grenseverdi for antatt konnektivitet mellom vernområdene", paste("basert på den funksjonale avstanden mellom dem i ", as.character(10^(nchar(as.integer(max(df_edgeremoval\$distance)))-2)), sep="")), ylab="Antall vernområder / nettverkskomponenter", lty=1, yaxt="n", yaxs="r", ylim=c(0, 100), yaxs="i", xaxs="i")
+#if (connectivity_cutoff>0) abline(v=connectivity_cutoff, col="red", lty=3)
+#ifelse(connectivity_cutoff>0, legend("topleft", c("Andel nettverkskomponenter", "Størrelsen av den største nettverkskomponenten", "(i andel populasjons størrelse)", "Antall forbindelser i nettverket", "Andel av diameter av nettverket", "Antatt grenseverdi for konnektivitet"), lty=c(1, 2, 0, 3, 4, 3), col=c("black", "black", "black", "black", "black", "red"), inset=0.005, bty="o", box.lty=0, bg="White"), legend("topleft", c("Andel nettverkskomponenter", "Størrelsen av den største nettverkskomponenten", "(i andel populasjons størrelse)", "Antall forbindelser i nettverket", "Andel av diameter av nettverket"), lty=c(1, 2, 0, 3, 4), col=c("black", "black", "black", "black", "black"), inset=0.005, bty="o", box.lty=0, bg="White"))
+#axis(2, seq.int(0, 100, 25), yaxs="r")
+#matplot(df_edgeremoval\$distance, df_edgeremoval\$cl_del_max_size*100/sum(V(g)\$pop_proxy), type="l", lty=2, yaxt="n", yaxs="i", add=TRUE)
+#lines(df_edgeremoval\$distance, dist_inv_ix*100/edges_n, type="l", lty=3)
+#lines(df_edgeremoval\$distance, df_edgeremoval\$cl_del_diam*100/diam_d, type="l", lty=4)
+##Lable axis 4!!!
+##axis(4, seq.int(0, maks_forbindelser_ud*maks_nodes/maks_forbindelser_ud, 100*maks_nodes/maks_forbindelser_ud), yaxs="i", labels=seq.int(0, maks_forbindelser_ud, 100), ylab="Antall forbindelser")
+#off <- dev.off(dev.cur())
+}
+
+########################################################################
+###Analysis on edge level
+########################################################################
+
+cat("Starting analysis on edge level")
+
+########################################################################
+######Calculate minimum spanning trees (MST) on the undirected graph with only direct edges
+###Minimum spanning tree weighted by maximum potential flow
+
+E(g_ud_d)\$mf_mst_ud <- as.integer(0)
+E(g_ud_d)[grep(TRUE, E(g_ud_d)\$con_id %in% E(minimum.spanning.tree(g_ud_d, weights=E(g_ud_d)\$mf_inv_u))\$con_id)]\$mf_mst_ud <- 1
+
+###Minimum spanning tree weighted by cost distance
+E(g_ud_d)\$cd_mst_ud <- as.integer(0)
+E(g_ud_d)[grep(TRUE, E(g_ud_d)\$con_id %in% E(minimum.spanning.tree(g_ud_d, weights=E(g_ud_d)\$cd_u))\$con_id)]\$cd_mst_ud <- 1
+
+E(g_ud_d)\$cf_mst_ud <- as.integer(0)
+E(g_ud_d)[grep(TRUE, E(g_ud_d)\$con_id %in% E(minimum.spanning.tree(g_ud_d, weights=E(g_ud_d)\$cf_inv_u))\$con_id)]\$cf_mst_ud <- 1
+
+######Calculate minimum spanning trees (MST) on the undirected graph with only direct edges shorter cost distance threshold
+###Minimum spanning tree weighted by maximum potential flow
+E(g_ud_d_cd)\$mf_mst_udc <- as.integer(0)
+E(g_ud_d_cd)[grep(TRUE, E(g_ud_d_cd)\$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)\$mf_inv_u))\$con_id)]\$mf_mst_udc <- 1
+
+
+###Minimum spanning tree weighted by cost distance
+E(g_ud_d_cd)\$cd_mst_udc <- as.integer(0)
+E(g_ud_d_cd)[grep(TRUE, E(g_ud_d_cd)\$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)\$cd_u))\$con_id)]\$cd_mst_udc <- 1
+
+###Minimum spanning tree weighted by competing potential flow
+E(g_ud_d_cd)\$cf_mst_udc <- as.integer(0)
+E(g_ud_d_cd)[grep(TRUE, E(g_ud_d_cd)\$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)\$cf_inv_u))\$con_id)]\$cf_mst_udc <- 1
+
+########################################################################
+###Identify bridges for the undirected graph
+E(g_ud)\$is_br_u <- as.integer(unlist(mclapply(1:(length(E(g_ud))), function(x) clusters(delete.edges(g_ud, E(g_ud)[x]))\$no-cl_no_d)))
+
+#Identify bridges for the undirected graph with only direct edges
+E(g_ud_d)\$is_br_ud <- as.integer(unlist(mclapply(1:(length(E(g_ud_d))), function(x) clusters(delete.edges(g_ud_d, E(g_ud_d)[x]))\$no-cl_no_d)))
+
+###Identify bridges for the undirected graph with only direct edges shorter cost distance threshold
+E(g_ud_d_cd)\$is_br_udc <- as.integer(unlist(mclapply(1:(length(E(g_ud_d_cd))), function(x) clusters(delete.edges(g_ud_d_cd, E(g_ud_d_cd)[x]))\$no-cl_no_d_cd)))
+
+########################################################################
+
+biconnected_components_d <- biconnected.components(g_ud)
+
+#Identify component edges (biconnected components) for the undirected graph with only direct edges
+E(g_ud)\$bc_e_u <- as.integer(0)
+for (c in 1:biconnected_components_d\$no) { 
+E(g_ud)\$bc_e_u[unlist(biconnected_components_d\$component_edges[c])] <- c
+}
+#Identify tree edges (biconnected components) for the undirected graph with only direct edges
+E(g_ud)\$bc_te_u <- as.integer(0)
+for (c in 1:biconnected_components_d\$no) {
+E(g_ud)\$bc_te_u[unlist(biconnected_components_d\$tree_edges[c])] <- c
+}
+
+biconnected_components_d <- biconnected.components(g_ud_d)
+#Identify component edges (biconnected components) for the undirected graph with only direct edges
+E(g_ud_d)\$bc_e_ud <- as.integer(0)
+for (c in 1:biconnected_components_d\$no) { 
+E(g_ud_d)\$bc_e_ud[unlist(biconnected_components_d\$component_edges[c])] <- c
+}
+#Identify tree edges (biconnected components) for the undirected graph with only direct edges
+E(g_ud_d)\$bc_te_ud <- as.integer(0)
+for (c in 1:biconnected_components_d\$no) {
+E(g_ud_d)\$bc_te_ud[unlist(biconnected_components_d\$tree_edges[c])] <- c
+}
+
+biconnected_components_d_cd <- biconnected.components(g_ud_d_cd)
+#Identify component edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold
+E(g_ud_d_cd)\$bc_e_udc <- as.integer(0)
+for (c in 1:biconnected_components_d_cd\$no) {
+E(g_ud_d_cd)\$bc_e_udc[unlist(biconnected_components_d_cd\$component_edges[c])] <- c
+}
+#Identify tree edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold
+E(g_ud_d_cd)\$bc_te_udc <- as.integer(0)
+for (c in 1:biconnected_components_d_cd\$no) {
+E(g_ud_d_cd)\$bc_te_udc[unlist(biconnected_components_d_cd\$tree_edges[c])] <- c
+}
+
+########################################################################
+###Calculate edge betweenness
+##Results of the internal edge.betweeness.estimate function are slightly different from the one used here (always larger values):
+##see: E(g_ud_d)\$cd_leb_ud <- edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=FALSE, lnbh_cutoff*connectivity_cutoff, weights=E(g_ud_d)\$cd_u)
+##But since the internal edge.betweeness.estimate function is not reasonable to use with other weights than cost distance, the workarounds are used nevertheless
+##Further investigation necessary!!!
+
+###Calculate edge betweenness on the entire undirected graph with only direct edges
+
+path_lengths <- shortest.paths(g_ud_d, v=V(g_ud_d), to=V(g_ud_d), weights=E(g_ud_d)\$cd_u)
+
+#weighted by cost distance
+E(g_ud_d)\$cd_eb_ud <- edge.betweenness(g_ud_d, e=E(g_ud_d), directed=FALSE, weights=E(g_ud_d)\$cd_u)
+E(g_ud_d)\$cd_leb_ud <- hist(unlist(mclapply(1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x,  V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)\$cd_u, output=c("epath")))), breaks=0:(length(E(g_ud_d))), plot=FALSE)\$counts/2
+#E(g_ud_d)\$cd_leb_ud <- edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=FALSE, lnbh_cutoff*connectivity_cutoff, weights=E(g_ud_d)\$cd_u)
+
+#weighted by maximum potential flow
+E(g_ud_d)\$mf_eb_ud <- edge.betweenness(g_ud_d, e=E(g_ud_d), directed=FALSE, weights=E(g_ud_d)\$mf_inv_u)
+E(g_ud_d)\$mf_leb_ud <- hist(unlist(mclapply(1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x,  V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)\$mf_inv_u, output=c("epath")))), breaks=0:(length(E(g_ud_d))), plot=FALSE)\$counts/2
+
+##weighted by competing potential flow
+E(g_ud_d)\$cf_eb_ud <- edge.betweenness(g_ud_d, e=E(g_ud_d), directed=FALSE, weights=E(g_ud_d)\$cf_inv_u)
+E(g_ud_d)\$cf_leb_ud <- hist(unlist(mclapply(1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x,  V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)\$cf_inv_u, output=c("epath")))), breaks=0:(length(E(g_ud_d))), plot=FALSE)\$counts/2
+
+########################################################################
+###Calculate edge betweenness on the entire undirected graph with only direct edges shorter cost distance threshold
+
+path_lengths_cd <- shortest.paths(g_ud_d_cd, v=V(g_ud_d_cd), to=V(g_ud_d_cd), weights=E(g_ud_d_cd)\$cd_u)
+
+#weighted by cost distance
+E(g_ud_d_cd)\$cd_eb_udc <- edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=FALSE, weights=E(g_ud_d_cd)\$cd_u)
+E(g_ud_d_cd)\$cd_leb_udc <- hist(unlist(mclapply(1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x,  V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)\$cd_u, output=c("epath")))), breaks=0:(length(E(g_ud_d_cd))), plot=FALSE)\$counts/2
+
+
+#weighted by maximum potential flow
+E(g_ud_d_cd)\$mf_eb_udc <- edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=FALSE, weights=E(g_ud_d_cd)\$mf_inv_u)
+E(g_ud_d_cd)\$mf_leb_udc <- hist(unlist(mclapply(1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x,  V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)\$mf_inv_u, output=c("epath")))), breaks=0:(length(E(g_ud_d_cd))), plot=FALSE)\$counts/2
+
+##weighted by competing potential flow
+E(g_ud_d_cd)\$cf_eb_udc <- edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=FALSE, weights=E(g_ud_d_cd)\$cf_inv_u)
+E(g_ud_d_cd)\$cf_leb_udc <- hist(unlist(mclapply(1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x,  V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)\$cf_inv_u, output=c("epath")))), breaks=0:(length(E(g_ud_d_cd))), plot=FALSE)\$counts/2
+
+#########################################################################
+#######Calculate edge betweenness community on the undirected graph
+##weighted by competing potential flow
+ebc <- edge.betweenness.community(g_ud, weights=E(g_ud)\$cf_inv_u, directed=FALSE, edge.betweenness=TRUE, merges=TRUE, bridges=TRUE, modularity=TRUE, membership=TRUE)
+
+E(g_ud)\$cf_iebc_v <- ebc\$edge.betweenness
+#cf_iebc_m <- ebc\$merges
+E(g_ud)\$cf_iebc_r <- ebc\$removed.edges
+E(g_ud)\$cf_iebc_b <- unlist(mclapply(1:length(E(g_ud)), function(x) ifelse(x %in% ebc\$bridges, 1, 0)))
+E(g_ud)\$cf_iebc_c <- crossing(ebc, g_ud)
+
+#V(g)\$cf_iebc_mo <- ebc\$modularity
+V(g)\$cf_iebc_me <- ebc\$membership
+
+com_struct <- as.character(cutat(ebc, cl_no_d))
+for (i in (cl_no_d+1):(cl_no_d+cl_thres)) {
+com_struct <- paste(com_struct, cutat(ebc, i), sep=';')
+}
+V(g)\$cf_iebc_cs <- as.character(com_struct)
+V(g)\$cf_iebc_cl <- cutat(ebc, (cl_no_d+cl_thres))
+
+cf_iebc_u_mod <- modularity(ebc)
+com_no_u <- length(ebc)
+com_sizes_u <- sizes(ebc)
+com_sizes_u_names <- paste("Size of comumity ", names(sizes(ebc)), " (at maximum modularity score (from iebc))", sep="")
+
+#The following workaround was written when igraph did not support weights when calculating edge.betweenness.community
+#The results of this procedure differ significantly from the (new) edge.betweenness.community in igraph 0.6-2 with weights (above)
+#In a first visual the results of the following procedure look more reasonable, but further investigation is necessary!!!
+
+n <- as.integer(1)
+ebc_con_id_u <- as.integer(0)
+ebc_rank <- as.integer(0)
+ebc_value <- as.integer(0)
+ebc_clust <- as.integer(0)
+V(g_ud)\$cf_ebc_cs <- as.character(clusters(g_ud)\$membership)
+del_eb_g <- g_ud
+
+for (edge in 1:length(E(del_eb_g))) {
+##Recalculate edge betweenness for the remaining edges
+E(del_eb_g)\$cf_eb_ud <- edge.betweenness(del_eb_g, e=E(del_eb_g), directed=FALSE, weights=E(del_eb_g)\$cf_inv_u)
+##Identify edge with highest edge betweenness value
+idx <- sort.int(E(del_eb_g)\$cf_eb_ud, decreasing=TRUE, na.last=NA, index.return=TRUE)\$ix[1]
+##Save edge betweenness value for this edge
+ebc_con_id_u[n] <- E(del_eb_g)\$con_id_u[idx]
+ebc_rank[n] <- n
+ebc_clust[n] <- clusters(del_eb_g)\$no
+if((ebc_clust[n]-cl_no_d)<(cl_no_d+cl_thres)) {if(n>1) {if(ebc_clust[n]>ebc_clust[n-1]) V(g_ud)\$cf_ebc_cs <- paste(V(g_ud)\$cf_ebc_cs, clusters(del_eb_g)\$membership, sep=";")}}
+ebc_value[n] <- E(del_eb_g)\$cf_eb_ud[idx]
+##Delete edge with highest edge betweenness value
+del_eb_g <- delete.edges(del_eb_g, E(del_eb_g)[idx])
+n <- n+1
+}
+
+V(g_ud)\$cf_ebc_cl <- unlist(mclapply(1:length(V(g_ud)), function(x) strsplit(V(g_ud)\$cf_ebc_cs, ";")[[x]][(cl_no_d+cl_thres)]))
+
+#ebc_df <- merge(data.frame(con_id=E(g)\$con_id, con_id_u=E(g)\$con_id_u), data.frame(con_id_u=ebc_con_id_u, ebc_rank=ebc_rank, ebc_clust=ebc_clust, ebc_value=ebc_value), all.x=TRUE)
+idx <- sort.int(ebc_con_id_u, index.return=TRUE)\$ix
+E(g_ud)\$cf_ebc_v <- ebc_value[idx]
+E(g_ud)\$cf_ebc_r <- ebc_rank[idx]
+E(g_ud)\$cf_ebc_c <- ebc_clust[idx]
+
+E(g_ud)\$cf_ebc_vi <- strftime(Sys.time()+E(g_ud)\$cf_ebc_r)
+########################################################################
+######Calculate cluster membership
+V(g_ud_d)\$cl_ud <- clusters(g_ud_d)\$membership
+V(g_ud_d_cd)\$cl_udc <- clusters(g_ud_d_cd)\$membership
+
+#########################################################################
+###Calculate potential cluster connectors (based on cost distance threshold)
+cl_membership <- data.frame(patch_id=V(g_ud_d_cd)\$patch_id, cl=V(g_ud_d_cd)\$cl_udc)
+cl_pc_pre <- merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)\$from_p, to_p=E(g)\$to_p), cl_membership, by.x="from_p", by.y="patch_id"),  cl_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)\$from_p, to_p=E(g)\$to_p), cl_membership, by.x="from_p", by.y="patch_id"),  cl_membership, by.x="to_p", by.y="patch_id")\$id) ,]
+E(g)\$cl_pc <- ifelse(cl_pc_pre\$cl.x==cl_pc_pre\$cl.y,0,1)
+
+#########################################################################
+###Calculate community connectors (based on community structure from ebc)
+com_membership <- data.frame(patch_id=V(g_ud_d_cd)\$patch_id, com=V(g_ud)\$cf_ebc_cl)
+com_pc_pre <- merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)\$from_p, to_p=E(g)\$to_p), com_membership, by.x="from_p", by.y="patch_id"),  com_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)\$from_p, to_p=E(g)\$to_p), com_membership, by.x="from_p", by.y="patch_id"),  com_membership, by.x="to_p", by.y="patch_id")\$id) ,]
+E(g)\$cf_ebc_cc <- ifelse(com_pc_pre\$com.x==com_pc_pre\$com.y,0,1)
+
+
+cat("Starting analysis on vertex level")
+
+########################################################################
+###Analysis on vertex level
+########################################################################
+
+########################################################################
+######Calculate degree centrality
+#On the directed graph
+#V(g)\$deg_dir <- as.integer(degree(g, v=V(g), mode=c("in")))
+
+#On the entire undirected graph
+V(g_ud)\$deg_u <- as.integer(degree(g_ud, v=V(g_ud)))
+
+#On the undirected graph with edges shorter cost distance threshold
+V(g_ud_cd)\$deg_uc <- as.integer(degree(g_ud_cd, v=V(g_ud_cd)))
+
+#On the undirected graph with only direct edges
+V(g_ud_d)\$deg_ud <- as.integer(degree(g_ud_d, v=V(g_ud_d)))
+
+#On the undirected graph with only direct edges shorter cost distance threshold
+V(g_ud_d_cd)\$deg_udc <- as.integer(degree(g_ud_d_cd, v=V(g_ud_d_cd)))
+
+########################################################################
+######Calculate closeness centrality
+V(g_ud_d)\$cd_cl_ud <- closeness(g_ud_d, vids=V(g_ud_d), weights = E(g_ud_d)\$cd_u, normalized = FALSE)
+V(g_ud_d)\$mf_cl_ud <- closeness(g_ud_d, vids=V(g_ud_d), weights = E(g_ud_d)\$mf_inv_u, normalized = FALSE)
+V(g_ud_d)\$cf_cl_ud <- closeness(g_ud_d, vids=V(g_ud_d), weights = E(g_ud_d)\$cf_inv_u, normalized = FALSE)
+
+V(g_ud_d_cd)\$cd_cl_udc <- closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)\$cd_u, normalized = FALSE)
+V(g_ud_d_cd)\$mf_cl_udc <- closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)\$mf_inv_u, normalized = FALSE)
+V(g_ud_d_cd)\$cf_cl_udc <- closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)\$cf_inv_u, normalized = FALSE)
+
+########################################################################
+######Calculate biconnected components
+
+##On the undirected graph with only direct edges
+#Number of new clusters when a vertex is deleted
+V(g_ud_d)\$art_ud <- as.integer(unlist(mclapply(1:(length(V(g_ud_d))), function(x) ifelse((clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))\$no-cl_no_d)<0,0,clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))\$no-cl_no_d))))
+#Vertex is articulation point
+V(g_ud_d)\$art_p_ud <- as.integer(0)
+V(g_ud_d)\$art_p_ud[biconnected_components_d\$articulation_points] <- as.integer(1)
+
+##On the undirected graph with only direct edges shorter cost distance threshold
+#Number of new clusters when a vertex is deleted
+V(g_ud_d_cd)\$art_udc <- unlist(mclapply(1:(length(V(g_ud_d_cd))), function(x) ifelse((clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))\$no-cl_no_d_cd)<0,0,clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))\$no-cl_no_d_cd)))
+#Vertex is articulation point
+V(g_ud_d_cd)\$art_p_udc <- as.integer(0)
+V(g_ud_d_cd)\$art_p_udc[biconnected_components_d_cd\$articulation_points] <- as.integer(1)
+
+########################################################################
+######Calculate eigenvector centrality
+#Calculate eigenvector for the entire directed graph
+
+#Calculate eigenvector centrality weighted by competing potential flow
+cf_evc_groups <- groupedData(cf ~ 1 | to_p, data=e, FUN=sum)
+cf_evc <- gsummary(cf_evc_groups, sum)
+df_cf_evc_d <- data.frame(patch_id=cf_evc\$to_p, cf_evc_d=cf_evc\$cf)
+
+#Calculate eigenvector centrality weighted by maximum potential flow
+mf_o_evc_groups <- groupedData(mf_o ~ 1 | to_p, data=e, FUN=sum)
+mf_o_evc <- gsummary(mf_o_evc_groups, sum)
+df_mf_i_evc_d <- data.frame(patch_id=mf_o_evc\$to_p, mf_evc_d=mf_o_evc\$mf_o)
+
+df_evc_d <- merge(df_cf_evc_d, df_mf_i_evc_d)
+
+######Calculate eigenvector centrality on the directed graph with edges shorter cost distance threshold
+#Calculate eigenvector centrality weighted by competing potential flow
+evc_udc_pre <- data.frame(to_p=e\$to_p[grep(TRUE, e\$cd_u<connectivity_cutoff)], mf_o=e\$mf_o[grep(TRUE, e\$cd_u<connectivity_cutoff)], cf=e\$cf[grep(TRUE, e\$cd_u<connectivity_cutoff)])
+
+cf_evc_udc_groups <- groupedData(cf ~ 1 | to_p, data=evc_udc_pre, FUN=sum)
+cf_evc_udc <- gsummary(cf_evc_udc_groups, sum)
+df_cf_evc_udc <- data.frame(patch_id=cf_evc_udc\$to_p, cf_evc=cf_evc_udc\$cf)
+
+#Calculate eigenvector centrality weighted by maximum potential flow
+mf_o_evc_udc_groups <- groupedData(mf_o ~ 1 | to_p, data=evc_udc_pre, FUN=sum)
+mf_o_evc_udc <- gsummary(mf_o_evc_udc_groups, sum)
+df_mf_o_evc_udc <- data.frame(patch_id=mf_o_evc_udc\$to_p, mf_evc=mf_o_evc_udc\$mf_o)
+
+df_evc_cd <- merge(df_cf_evc_udc, df_mf_o_evc_udc)
+
+###Add eigenvector centrality as a vertex attribute
+for (p in 1:length(V(g))) {
+	V(g)\$mf_evc_d[p] <- ifelse(length(grep(TRUE, as.integer(as.character(df_evc_d\$patch_id))==as.integer(V(g)\$patch_id[p])))<1,0,df_evc_d\$mf_evc[grep(TRUE, as.integer(as.character(df_evc_d\$patch_id))==as.integer(V(g)\$patch_id[p]))])
+	V(g)\$cf_evc_d[p] <- ifelse(length(grep(TRUE, as.integer(as.character(df_evc_d\$patch_id))==as.integer(V(g)\$patch_id[p])))<1,0,df_evc_d\$cf_evc[grep(TRUE, as.integer(as.character(df_evc_d\$patch_id))==as.integer(V(g)\$patch_id[p]))])
+	V(g)\$mf_evc_cd[p] <- ifelse(length(grep(TRUE, as.integer(as.character(df_evc_cd\$patch_id))==as.integer(V(g)\$patch_id[p])))<1,0,df_evc_cd\$mf_evc[grep(TRUE, as.integer(as.character(df_evc_cd\$patch_id))==as.integer(V(g)\$patch_id[p]))])
+	V(g)\$cf_evc_cd[p] <- ifelse(length(grep(TRUE, as.integer(as.character(df_evc_cd\$patch_id))==as.integer(V(g)\$patch_id[p])))<1,0,df_evc_cd\$cf_evc[grep(TRUE, as.integer(as.character(df_evc_cd\$patch_id))==as.integer(V(g)\$patch_id[p]))])
+}
+
+
+########################################################################
+###Calculate vertex betweenness
+##Results of the internal betweeness.estimate function are slightly different from the one used here (always larger values)
+##see: V(g_ud_d)\$cd_lvb_ud <- betweenness.estimate(g_ud_d, vids = V(g_ud_d), directed = FALSE, lnbh_cutoff*connectivity_cutoff, weights = E(g_ud_d)\$cd_u)
+##But since the internal betweeness.estimate function is not reasonable to use with other weights than cost distance, the workarounds are used nevertheless
+##Further investigation necessary!!!
+
+###Calculate vertex betweenness on the undirected graph with only direct edges
+V(g_ud_d)\$cd_vb_ud <- betweenness(g_ud_d, v=V(g_ud_d), directed = FALSE, weights = E(g_ud_d)\$cd_u)
+vsp_cd_ud <- function(x) unlist(get.shortest.paths(g_ud_d, x,  V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)\$cd_u, output=c("vpath")))
+V(g_ud_d)\$cd_lvb_ud <- hist(unlist(mclapply(1:length(V(g_ud_d)), function(x) vsp_cd_ud(x)[grep(FALSE, 1:length(vsp_cd_ud(x)) %in% append(c(1, length(vsp_cd_ud(x))), append(grep(TRUE, vsp_cd_ud(x)[1:length(vsp_cd_ud(x))]==x), grep(TRUE, vsp_cd_ud(x)[1:length(vsp_cd_ud(x))]==x)-1)))])), breaks=0:length(V(g_ud_d)), plot=FALSE)\$counts/2
+
+##weighted by maximum potential flow
+V(g_ud_d)\$mf_vb_ud <- betweenness(g_ud_d, v=V(g_ud_d), directed = TRUE, weights = E(g_ud_d)\$mf_inv_u)
+vsp_mf_u <- function(x) unlist(get.shortest.paths(g_ud_d, x,  V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)\$mf_inv_u, output=c("vpath")))
+V(g_ud_d)\$mf_lvb_ud <- hist(unlist(mclapply(1:length(V(g_ud_d)), function(x) vsp_mf_u(x)[grep(FALSE, 1:length(vsp_mf_u(x)) %in% append(c(1, length(vsp_mf_u(x))), append(grep(TRUE, vsp_mf_u(x)[1:length(vsp_mf_u(x))]==x), grep(TRUE, vsp_mf_u(x)[1:length(vsp_mf_u(x))]==x)-1)))])), breaks=0:length(V(g_ud_d)), plot=FALSE)\$counts/2
+
+##weighted by competing potential flow
+V(g_ud_d)\$cf_vb_ud <- betweenness(g_ud_d, v=V(g_ud_d), directed = TRUE, weights = E(g_ud_d)\$cf_inv_u)
+vsp_cf_u <- function(x) unlist(get.shortest.paths(g_ud_d, x,  V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)\$cf_inv_u, output=c("vpath")))
+V(g_ud_d)\$cf_lvb_ud <- hist(unlist(mclapply(1:length(V(g_ud_d)), function(x) vsp_cf_u(x)[grep(FALSE, 1:length(vsp_cf_u(x)) %in% append(c(1, length(vsp_cf_u(x))), append(grep(TRUE, vsp_cf_u(x)[1:length(vsp_cf_u(x))]==x), grep(TRUE, vsp_cf_u(x)[1:length(vsp_cf_u(x))]==x)-1)))])), breaks=0:length(V(g_ud_d)), plot=FALSE)\$counts/2
+
+###Calculate vertex betweenness on the undirected graph with only direct edges shorter connectivity cutoff
+V(g_ud_d_cd)\$cd_vb_udc <- betweenness(g_ud_d_cd, v=V(g_ud_d_cd), directed = FALSE, weights = E(g_ud_d_cd)\$cd_u)
+vsp_cd_udc <- function(x) unlist(get.shortest.paths(g_ud_d_cd, x,  V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)\$cd_u, output=c("vpath")))
+V(g_ud_d_cd)\$cd_lvb_udc <- hist(unlist(mclapply(1:length(V(g_ud_d_cd)), function(x) vsp_cd_udc(x)[grep(FALSE, 1:length(vsp_cd_udc(x)) %in% append(c(1, length(vsp_cd_udc(x))), append(grep(TRUE, vsp_cd_udc(x)[1:length(vsp_cd_udc(x))]==x), grep(TRUE, vsp_cd_udc(x)[1:length(vsp_cd_udc(x))]==x)-1)))])), breaks=0:length(V(g_ud_d_cd)), plot=FALSE)\$counts/2
+
+##weighted by maximum potential flow
+V(g_ud_d_cd)\$mf_vb_udc <- betweenness(g_ud_d_cd, v=V(g_ud_d_cd), directed = TRUE, weights = E(g_ud_d_cd)\$mf_i_inv_ud)
+vsp_mf_uc <- function(x) unlist(get.shortest.paths(g_ud_d_cd, x,  V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)\$mf_inv_u, output=c("vpath")))
+V(g_ud_d_cd)\$mf_lvb_udc <- hist(unlist(mclapply(1:length(V(g_ud_d_cd)), function(x) vsp_mf_uc(x)[grep(FALSE, 1:length(vsp_mf_uc(x)) %in% append(c(1, length(vsp_mf_uc(x))), append(grep(TRUE, vsp_mf_uc(x)[1:length(vsp_mf_uc(x))]==x), grep(TRUE, vsp_mf_uc(x)[1:length(vsp_mf_uc(x))]==x)-1)))])), breaks=0:length(V(g_ud_d_cd)), plot=FALSE)\$counts/2
+
+##weighted by competing potential flow
+V(g_ud_d_cd)\$cf_vb_udc <- betweenness(g_ud_d_cd, v=V(g_ud_d_cd), directed = TRUE, weights = E(g_ud_d_cd)\$cf_inv_u)
+vsp_cf_uc <- function(x) unlist(get.shortest.paths(g_ud_d_cd, x,  V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)\$cf_inv_u, output=c("vpath")))
+V(g_ud_d_cd)\$cf_lvb_udc <- hist(unlist(mclapply(1:length(V(g_ud_d_cd)), function(x) vsp_cf_uc(x)[grep(FALSE, 1:length(vsp_cf_uc(x)) %in% append(c(1, length(vsp_cf_uc(x))), append(grep(TRUE, vsp_cf_uc(x)[1:length(vsp_cf_uc(x))]==x), grep(TRUE, vsp_cf_uc(x)[1:length(vsp_cf_uc(x))]==x)-1)))])), breaks=0:length(V(g_ud_d_cd)), plot=FALSE)\$counts/2
+
+########################################################################
+###Calculate neighborhood size
+V(g)\$nbh_s_uc <- as.integer(neighborhood.size(g_ud_cd, 1, nodes=V(g), mode=c("all")))
+###Calculate local neighborhood size
+V(g)\$nbh_sl_uc <- as.integer(neighborhood.size(g_ud_cd, lnbh_cutoff, nodes=V(g), mode=c("all")))
+
+##############################################
+#Export network overview measures
+network_overview_measures <- data.frame(Measure=c("Measure", "Number of vertices", "Number of edges (undirected)", "Number of direct edges (undirected)", "Number of edges shorter than cost distance threshold (undirected)", "Number of direct edges shorter than cost distance threshold (undirected)", "Number of clusters of the entire graph", "Number of clusters of the graph with only edges shorter cost distance threshold", "Size of the largest cluster of the entire graph", "Size of the largest cluster of the graph with only edges shorter cost distance threshold", "Average size of the clusters of the entire graph", "Average size of the clusters of the graph with only edges shorter cost distance threshold", "Diameter of the entire graph (undirected)", "Diameter of the graph with only direct edges (undirected)", "Diameter of the graph with only edges shorter cost distance threshold", "Density of the entire graph (directed)", "Density of the entire graph (undirected)", "Density of 
 the graph with only direct edges (undirected)", "Density of the graph with only edges shorter cost distance threshold", "Modularity (from iebc) of the entire graph (undirected) weighted by cf", "Number of communities (at maximum modularity score (from iebc)) of the entire (undirected) graph weighted by cf", com_sizes_u_names), Value=c("Value", vertices_n, edges_n, edges_n_d, edges_n_cd, edges_n_d_cd, cl_no_d, cl_no_d_cd, cl_max_size_d, cl_max_size_d_cd, cl_mean_size_d, cl_mean_size_d_cd, diam, diam_d, diam_d_cd, density, density_u, density_ud, density_udc, cf_iebc_u_mod, com_no_u, com_sizes_u))
+con_network_export <- file(network_output, open="wt")
+write.csv(network_overview_measures, con_network_export, row.names=FALSE, quote=FALSE)
+close(con_network_export)
+
+
+###############################################################
+#####Merge vertexmeasures in a dataframe and save it
+
+###Create initial export data frame
+
+vertex_export_df_ud <- as.data.frame(1:length(V(g_ud)))
+
+####Adjust vertex-attributes first
+
+g_ud <- remove.vertex.attribute(g_ud, "centroid_x")
+g_ud <- remove.vertex.attribute(g_ud, "centroid_y")
+g_ud <- remove.vertex.attribute(g_ud, "pop_proxy")
+
+vertex_attribute_list <- list.vertex.attributes(g_ud)
+if(length(vertex_attribute_list)>0) {
+	for (vl in vertex_attribute_list) {
+		vertex_export_df_ud <- as.data.frame(cbind(vertex_export_df_ud, get.vertex.attribute(g_ud, vl)))
+	}
+
+	vertex_export_df_ud <- vertex_export_df_ud[2:length(vertex_export_df_ud)]
+	names(vertex_export_df_ud) <- vertex_attribute_list
+}
+
+###Create initial export data frame
+
+vertex_export_df_ud_d <- as.data.frame(1:length(V(g_ud_d)))
+
+####Adjust vertex-attributes first
+
+g_ud_d <- remove.vertex.attribute(g_ud_d, "centroid_x")
+g_ud_d <- remove.vertex.attribute(g_ud_d, "centroid_y")
+g_ud_d <- remove.vertex.attribute(g_ud_d, "pop_proxy")
+
+vertex_attribute_list <- list.vertex.attributes(g_ud_d)
+if(length(vertex_attribute_list)>0) {
+	for (vl in vertex_attribute_list) {
+		vertex_export_df_ud_d <- as.data.frame(cbind(vertex_export_df_ud_d, get.vertex.attribute(g_ud_d, vl)))
+	}
+	vertex_export_df_ud_d <- vertex_export_df_ud_d[2:length(vertex_export_df_ud_d)]
+	names(vertex_export_df_ud_d) <- vertex_attribute_list
+}
+
+###Create initial export data frame
+
+vertex_export_df_ud_d_cd <- as.data.frame(1:length(V(g_ud_d_cd)))
+
+####Adjust vertex-attributes first
+
+g_ud_d_cd <- remove.vertex.attribute(g_ud_d_cd, "centroid_x")
+g_ud_d_cd <- remove.vertex.attribute(g_ud_d_cd, "centroid_y")
+g_ud_d_cd <- remove.vertex.attribute(g_ud_d_cd, "pop_proxy")
+
+vertex_attribute_list <- list.vertex.attributes(g_ud_d_cd)
+
+if(length(vertex_attribute_list)>0) {
+	for (vl in vertex_attribute_list) {
+		vertex_export_df_ud_d_cd <- as.data.frame(cbind(vertex_export_df_ud_d_cd, get.vertex.attribute(g_ud_d_cd, vl)))
+	}
+	vertex_export_df_ud_d_cd <- vertex_export_df_ud_d_cd[2:length(vertex_export_df_ud_d_cd)]
+	names(vertex_export_df_ud_d_cd) <- vertex_attribute_list
+}
+###Create initial export data frame
+
+vertex_export_df_d <- as.data.frame(1:length(V(g)))
+
+
+V(g)\$WKT <- paste('POINT(', V(g)\$centroid_x, " ", V(g)\$centroid_y, ')', sep="")
+
+g <- remove.vertex.attribute(g, "centroid_x")
+g <- remove.vertex.attribute(g, "centroid_y")
+
+vertex_attribute_list <- list.vertex.attributes(g)
+if(length(vertex_attribute_list)>0) {
+
+	for (vl in vertex_attribute_list) {
+		vertex_export_df_d <- as.data.frame(cbind(vertex_export_df_d, get.vertex.attribute(g, vl)))
+	}
+	vertex_export_df_d <- vertex_export_df_d[2:length(vertex_export_df_d)]
+	names(vertex_export_df_d) <- vertex_attribute_list
+}
+
+vertex_export_df <- merge(merge(merge(vertex_export_df_d, vertex_export_df_ud, by="patch_id"), vertex_export_df_ud_d, by="patch_id"), vertex_export_df_ud_d_cd, by="patch_id")
+#df_1 <- data.frame(patch_id=V(g)\$patch_id, WKT=paste('POINT(', V(g)\$centroid_x, " ", V(g)\$centroid_y, ')', sep=""), pop_proxy=V(g)\$pop_proxy, cost_distance_vb=V(g)\$cost_distance_vb, cost_distance_vb_cutoff=V(g)\$cost_distance_vb_cutoff, mf_vb=V(g)\$mf_vb, mf_vb_cutoff=V(g)\$mf_vb_cutoff, cf_vb=V(g)\$cf_vb, cf_vb_cutoff=V(g)\$cf_vb_cutoff)
+#df_2 <- data.frame(patch_id=V(g)\$patch_id, cost_distance_vb_ud=V(g)\$cost_distance_vb_ud, cost_distance_vb_cutoff_ud=V(g)\$cost_distance_vb_cutoff_ud, mf_vb_ud=V(g)\$mf_vb_ud, mf_vb_cutoff_ud=V(g)\$mf_vb_cutoff_ud, cf_vb_ud=V(g)\$cf_vb_ud, cf_vb_cutoff_ud=V(g)\$cf_vb_cutoff_ud)
+#vertex_export_df <- merge(merge(df_1, df_2), df_evc)
+con_vertex_export <- file(vertex_output, open="wt")
+write.csv(vertex_export_df, con_vertex_export, row.names=FALSE)
+close(con_vertex_export)
+
+####Prepare VRT-files according to requested measures for vertex output
+vertex_vrt_var <- as.character()
+
+vertex_attribute_list <- sort(names(vertex_export_df))
+
+int_width <- 0
+is_int <- 0
+c_width <- 0
+for (vl in vertex_attribute_list) {
+	attrib <- na.omit(vertex_export_df[[match(vl, names(vertex_export_df))]])
+	if(is.factor(attrib)) {c_width <- max(nchar(as.character(attrib)))}
+	if(storage.mode(attrib)=="character") {c_width <- max(nchar(attrib))}
+	if(is.factor(attrib)==FALSE && storage.mode(attrib)!="character" && storage.mode(attrib)!="logical") {is_int <- ifelse(max(nchar(attrib)-nchar(round(attrib)))==0, 1, 0)}
+	if(is_int==1) {int_width <- max(nchar(round(attrib)))}
+	if(storage.mode(attrib)=="logical") {is_int <- 1 ; int_width <- 1}
+	vertex_vrt_var <- append(vertex_vrt_var, paste("<Field name=\"", vl, "\" src=\"", vl, "\" type=\"", ifelse(c_width>0, paste("String\" width=\"", c_width, sep=''), ifelse(is_int==0, "Real", paste("Integer\" width=\"", int_width, sep=''))), "\"/>", sep=''))
+	int_width <- 0
+	is_int <- 0
+	c_width <- 0
+}
+
+vertex_vrt <- c("<OGRVRTDataSource>", '<OGRVRTLayer name="vertex_measures">', paste("<SrcDataSource>", vertex_output, "</SrcDataSource>", sep=''), "<SrcLayer>vertex_measures</SrcLayer>", paste("<SrcSQL>SELECT", toString(names(vertex_export_df)), "FROM vertex_measures</SrcSQL>", sep=' '))
+vertex_vrt <- append(vertex_vrt, vertex_vrt_var)
+vertex_vrt <- append(vertex_vrt, c("</OGRVRTLayer>", "</OGRVRTDataSource>"))
+
+con_vertex_vrt <- file(vertex_output_vrt, open="wt")
+write.table(vertex_vrt, con_vertex_vrt, append = FALSE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)
+close(con_vertex_vrt)
+
+#########################################################################
+
+###Create initial export data frame for the undirected graph
+
+edge_export_df_ud <- as.data.frame(E(g_ud)\$con_id_u)
+
+####Adjust edge attributes for the undirected graphs first
+g_ud <- remove.edge.attribute(g_ud, "con_id")
+g_ud <- remove.edge.attribute(g_ud, "con_id_u")
+g_ud <- remove.edge.attribute(g_ud, "from_p")
+g_ud <- remove.edge.attribute(g_ud, "from_pop")
+g_ud <- remove.edge.attribute(g_ud, "to_p")
+g_ud <- remove.edge.attribute(g_ud, "to_pop")
+g_ud <- remove.edge.attribute(g_ud, "from_p_centroid_x")
+g_ud <- remove.edge.attribute(g_ud, "from_p_centroid_y")
+g_ud <- remove.edge.attribute(g_ud, "to_p_centroid_x")
+g_ud <- remove.edge.attribute(g_ud, "to_p_centroid_y")
+g_ud <- remove.edge.attribute(g_ud, "cost_distance")
+g_ud <- remove.edge.attribute(g_ud, "cd_u")
+g_ud <- remove.edge.attribute(g_ud, "distance_weight_e")
+g_ud <- remove.edge.attribute(g_ud, "distance_weight_e_ud")
+g_ud <- remove.edge.attribute(g_ud, "mf_o")
+g_ud <- remove.edge.attribute(g_ud, "mf_i")
+g_ud <- remove.edge.attribute(g_ud, "mf_o_inv")
+g_ud <- remove.edge.attribute(g_ud, "mf_i_inv")
+g_ud <- remove.edge.attribute(g_ud, "mf_u")
+g_ud <- remove.edge.attribute(g_ud, "mf_inv_u")
+g_ud <- remove.edge.attribute(g_ud, "cf")
+g_ud <- remove.edge.attribute(g_ud, "cf_inv")
+g_ud <- remove.edge.attribute(g_ud, "cf_u")
+g_ud <- remove.edge.attribute(g_ud, "cf_inv_u")
+
+edge_attribute_list <- list.edge.attributes(g_ud)
+
+if(length(edge_attribute_list )>0) {
+	for (el in edge_attribute_list) {
+		edge_export_df_ud <- as.data.frame(cbind(edge_export_df_ud, get.edge.attribute(g_ud, el)))
+	}
+	names(edge_export_df_ud) <- append("con_id_u", edge_attribute_list)
+}
+
+###Create initial export data frame for the undirected graph with only direct edges
+
+edge_export_df_ud_d <- as.data.frame(E(g_ud_d)\$con_id_u)
+
+####Adjust edge attributes for the undirected graph with only direct edges
+g_ud_d <- remove.edge.attribute(g_ud_d, "con_id")
+g_ud_d <- remove.edge.attribute(g_ud_d, "con_id_u")
+g_ud_d <- remove.edge.attribute(g_ud_d, "from_p")
+g_ud_d <- remove.edge.attribute(g_ud_d, "from_pop")
+g_ud_d <- remove.edge.attribute(g_ud_d, "to_p")
+g_ud_d <- remove.edge.attribute(g_ud_d, "to_pop")
+g_ud_d <- remove.edge.attribute(g_ud_d, "from_p_centroid_x")
+g_ud_d <- remove.edge.attribute(g_ud_d, "from_p_centroid_y")
+g_ud_d <- remove.edge.attribute(g_ud_d, "to_p_centroid_x")
+g_ud_d <- remove.edge.attribute(g_ud_d, "to_p_centroid_y")
+g_ud_d <- remove.edge.attribute(g_ud_d, "cost_distance")
+g_ud_d <- remove.edge.attribute(g_ud_d, "cd_u")
+g_ud_d <- remove.edge.attribute(g_ud_d, "distance_weight_e")
+g_ud_d <- remove.edge.attribute(g_ud_d, "distance_weight_e_ud")
+g_ud_d <- remove.edge.attribute(g_ud_d, "mf_o")
+g_ud_d <- remove.edge.attribute(g_ud_d, "mf_i")
+g_ud_d <- remove.edge.attribute(g_ud_d, "mf_o_inv")
+g_ud_d <- remove.edge.attribute(g_ud_d, "mf_i_inv")
+g_ud_d <- remove.edge.attribute(g_ud_d, "mf_u")
+g_ud_d <- remove.edge.attribute(g_ud_d, "mf_inv_u")
+g_ud_d <- remove.edge.attribute(g_ud_d, "cf")
+g_ud_d <- remove.edge.attribute(g_ud_d, "cf_inv")
+g_ud_d <- remove.edge.attribute(g_ud_d, "cf_u")
+g_ud_d <- remove.edge.attribute(g_ud_d, "cf_inv_u")
+g_ud_d <- remove.edge.attribute(g_ud_d, "isshort")
+g_ud_d <- remove.edge.attribute(g_ud_d, "isshort_cd")
+g_ud_d <- remove.edge.attribute(g_ud_d, "isshort_mf")
+g_ud_d <- remove.edge.attribute(g_ud_d, "isshort_cf")
+
+edge_attribute_list <- list.edge.attributes(g_ud_d)
+
+if(length(edge_attribute_list )>0) {
+	for (el in edge_attribute_list) {
+		edge_export_df_ud_d <- as.data.frame(cbind(edge_export_df_ud_d, get.edge.attribute(g_ud_d, el)))
+	}
+	names(edge_export_df_ud_d) <- append("con_id_u", edge_attribute_list)
+}
+###Create initial export data frame for the undirected graph with only direct edges
+
+edge_export_df_ud_d_cd <- as.data.frame(E(g_ud_d_cd)\$con_id_u)
+
+####Adjust edge attributes for the undirected graph with only direct edges
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "con_id")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "con_id_u")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "from_p")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "from_pop")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "to_p")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "to_pop")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "from_p_centroid_x")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "from_p_centroid_y")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "to_p_centroid_x")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "to_p_centroid_y")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cost_distance")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cd_u")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "distance_weight_e")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "distance_weight_e_ud")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_o")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_i")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_o_inv")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_i_inv")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_u")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_inv_u")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf_inv")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf_u")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf_inv_u")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort_cd")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort_mf")
+g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort_cf")
+
+edge_attribute_list <- list.edge.attributes(g_ud_d_cd)
+
+if(length(edge_attribute_list )>0) {
+	for (el in edge_attribute_list) {
+		edge_export_df_ud_d_cd <- as.data.frame(cbind(edge_export_df_ud_d_cd, get.edge.attribute(g_ud_d_cd, el)))
+	}
+	names(edge_export_df_ud_d_cd) <- append("con_id_u", edge_attribute_list)
+}
+
+
+###Create initial export data frame for the directed graph
+edge_export_df <- as.data.frame(1:length(E(g)))
+
+####Adjust edge attributes for the directed graph first
+
+E(g)\$WKT <- as.character(paste('LINESTRING(', E(g)\$from_p_centroid_x, " ", E(g)\$from_p_centroid_y, ", ", E(g)\$to_p_centroid_x, " ", E(g)\$to_p_centroid_y, ')', sep=""))
+
+g <- remove.edge.attribute(g, "from_p_centroid_x")
+g <- remove.edge.attribute(g, "from_p_centroid_y")
+g <- remove.edge.attribute(g, "to_p_centroid_x")
+g <- remove.edge.attribute(g, "to_p_centroid_y")
+E(g)\$cd <- E(g)\$cost_distance
+E(g)\$distk <- E(g)\$distance_weight_e
+E(g)\$distk_u <- E(g)\$distance_weight_e_ud
+g <- remove.edge.attribute(g, "distance_weight_e")
+g <- remove.edge.attribute(g, "distance_weight_e_ud")
+g <- remove.edge.attribute(g, "cost_distance")
+
+edge_attribute_list <- list.edge.attributes(g)
+
+if(length(edge_attribute_list )>0) {
+	for (el in edge_attribute_list) {
+		edge_export_df <- as.data.frame(cbind(edge_export_df, get.edge.attribute(g, el)))
+	}
+	names(edge_export_df) <- append("id", edge_attribute_list)
+}
+
+export_df_list <-c("edge_export_df_ud", "edge_export_df_ud_d", "edge_export_df_ud_d_cd")
+for (df in export_df_list) {
+	if(length(names(get(df)))) {
+		edge_export_df_final <- merge(edge_export_df, get(df), all.x=TRUE, by="con_id_u", suffixes=c("_x", "_y"))
+		edge_export_df <- edge_export_df_final
+	}
+}
+
+#edge_export_df <- merge(merge(merge(merge(merge(e, data.frame(con_id=E(g)\$con_id, WKT=paste('LINESTRING(', E(g)\$from_p_centroid_x, " ", E(g)\$from_p_centroid_y, ", ", E(g)\$to_p_centroid_x, " ", E(g)\$to_p_centroid_y, ')', sep=""), is_direct=E(g)\$isshort, cost_distance_eb=E(g)\$cost_distance_eb, cost_distance_eb_cutoff=E(g)\$cost_distance_eb_cutoff, mf_eb=E(g)\$mf_eb, mf_eb_cutoff=E(g)\$mf_eb_cutoff, cf_eb=E(g)\$cf_eb, cf_eb_cutoff=E(g)\$cf_eb_cutoff), all.x=TRUE), cost_distance_eb_ud, all.x=TRUE), mf_eb_ud, all.x=TRUE), cf_eb_ud, all.x=TRUE), mst_df, all.x=TRUE)
+#export_df <- merge(merge(merge(merge(data.frame(con_id=E(g)\$con_id, con_id_u=E(g)\$con_id_u, centroid_connection=paste('LINESTRING(', E(g)\$from_p_centroid_x, " ", E(g)\$from_p_centroid_y, ", ", E(g)\$to_p_centroid_x, " ", E(g)\$to_p_centroid_y, ')', sep=""), cost_distance_eb=E(g)\$cost_distance_eb, cost_distance_eb_cutoff=E(g)\$cost_distance_eb_cutoff, mf_eb=E(g)\$mf_eb, mf_eb_cutoff=E(g)\$mf_eb_cutoff, cf_eb=E(g)\$cf_eb, cf_eb_cutoff=E(g)\$cf_eb_cutoff), cost_distance_eb, all.x=TRUE), mf_eb, all.x=TRUE), cf_eb, all.x=TRUE), mst_df, all.x=TRUE)
+#edge_export_df <- merge(merge(merge(merge(merge(e, data.frame(con_id=E(g)\$con_id, centroid_connection=paste('"LINESTRING(', E(g)\$from_p_centroid_x, " ", E(g)\$from_p_centroid_y, ", ", E(g)\$to_p_centroid_x, " ", E(g)\$to_p_centroid_y, ')"', sep="")), all.x=TRUE), cost_distance_eb, all.x=TRUE), mf_eb, all.x=TRUE), cf_eb, all.x=TRUE), mst_df, all.x=TRUE)
+con_edge_export <- file(edge_output, open="wt")
+write.csv(edge_export_df_final, con_edge_export, row.names=FALSE, na="NULL")
+close(con_edge_export)
+
+####Prepare VRT-files according to requested measures for edge output
+edge_vrt_var <- as.character()
+
+edge_attribute_list <- sort(names(edge_export_df_final))
+
+int_width <- 0
+is_int <- 0
+c_width <- 0
+for (el in edge_attribute_list) {
+	attrib <- na.omit(edge_export_df_final[[match(el, names(edge_export_df_final))]])
+	if(is.factor(attrib)) {c_width <- max(nchar(as.character(attrib)))}
+	if(storage.mode(attrib)=="character") {c_width <- max(nchar(attrib))}
+	if(is.factor(attrib)==FALSE && storage.mode(attrib)!="character" && storage.mode(attrib)!="logical") {is_int <- ifelse(max(nchar(attrib)-nchar(round(attrib)))==0, 1, 0)}
+	if(is_int==1) {int_width <- max(nchar(round(attrib)))}
+	if(storage.mode(attrib)=="logical") {is_int <- 1 ; int_width <- 1}
+	edge_vrt_var <- append(edge_vrt_var, paste("<Field name=\"", el, "\" src=\"", el, "\" type=\"", ifelse(c_width>0, paste("String\" width=\"", c_width, sep=''), ifelse(is_int==0, "Real", paste("Integer\" width=\"", int_width, sep=''))), "\"/>", sep=''))
+	int_width <- 0
+	is_int <- 0
+	c_width <- 0
+}
+
+edge_vrt <- c("<OGRVRTDataSource>", '<OGRVRTLayer name="edge_measures">', paste("<SrcDataSource>", edge_output, "</SrcDataSource>", sep=''), "<SrcLayer>edge_measures</SrcLayer>", paste("<SrcSQL>SELECT", toString(names(edge_export_df)), "FROM edge_measures</SrcSQL>", sep=' '))
+edge_vrt <- append(edge_vrt, edge_vrt_var)
+edge_vrt <- append(edge_vrt, c("</OGRVRTLayer>", "</OGRVRTDataSource>"))
+
+con_edge_vrt <- file(edge_output_vrt, open="wt")
+write.table(edge_vrt, con_edge_vrt, append = FALSE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)
+close(con_edge_vrt)
+
+#########################
+#CREATE .qml-files for edge measures visualistion in QGIS
+
+no_quantile <- 5
+colortable <- c("          <prop k=\"color\" v=\"215,25,28,255\"/>", "          <prop k=\"color\" v=\"253,174,97,255\"/>", "          <prop k=\"color\" v=\"255,255,191,255\"/>", "          <prop k=\"color\" v=\"166,217,106,255\"/>", "          <prop k=\"color\" v=\"26,150,65,255\"/>")
+colortable_bin <- c("          <prop k=\"color\" v=\"0,0,0,255\"/>")
+for (attribute in names(edge_export_df_final)) {
+ 
+#Skip id and geom columns
+if(attribute %in% c("id", "con_id", "con_id_u", "from_p", "to_p", "WKT", "cf_ebc_vi")) { next }
+
+st_mod <- storage.mode(unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))]))
+att_val <- unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))])
+
+#Write header
+qml <- c("<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>")
+qml <- append(qml, '<qgis version="1.8" minimumScale="0" maximumScale="1e+08" hasScaleBasedVisibilityFlag="0">')
+qml <- append(qml, '  <transparencyLevelInt>255</transparencyLevelInt>')
+
+if((max(att_val, na.rm=TRUE)-min(att_val, na.rm=TRUE)==1)) {
+
+#More header
+qml <- append(qml, paste('  <renderer-v2 attr="', attribute, '" symbollevels="0" type="categorizedSymbol">', sep=''))
+
+#Categories
+qml <- append(qml, "    <categories>")
+qml <- append(qml, "      <category symbol=\"0\" value=\"1\" label=\"\"/>")
+qml <- append(qml, "    </categories>")
+
+#Write symbols
+qml <- append(qml, "    <symbols>")
+qml <- append(qml, '      <symbol outputUnit="MM" alpha="1" type="line" name="0">')
+qml <- append(qml, "        <layer pass=\"0\" class=\"SimpleLine\" locked=\"0\">")
+qml <- append(qml, "          <prop k=\"capstyle\" v=\"square\"/>")
+qml <- append(qml, colortable_bin)
+qml <- append(qml, "          <prop k=\"customdash\" v=\"5;2\"/>")
+qml <- append(qml, "          <prop k=\"joinstyle\" v=\"bevel\"/>")
+qml <- append(qml, "          <prop k=\"offset\" v=\"0\"/>")
+qml <- append(qml, "          <prop k=\"penstyle\" v=\"solid\"/>")
+qml <- append(qml, "          <prop k=\"use_custom_dash\" v=\"0\"/>")
+qml <- append(qml, "          <prop k=\"width\" v=\"0.26\"/>")
+qml <- append(qml, "        </layer>")
+qml <- append(qml, "      </symbol>")
+}
+else {
+
+attribute_quantile <- quantile(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))], probs=seq(0, 1, 1/5), na.rm=TRUE, type=8)
+
+#Write more header
+qml <- append(qml, paste('  <renderer-v2 attr="', attribute, '" symbollevels="1" type="graduatedSymbol">', sep=''))
+
+#Write ranges
+ranges <- character()
+qml <- append(qml, "    <ranges>")
+for (quant in 1:no_quantile) {
+ranges <- append(ranges, paste("      <range symbol=\"", (quant-1), "\" lower=\"", attribute_quantile[quant], "\" upper=\"", attribute_quantile[(quant+1)], "\" label=\"", round(attribute_quantile[quant], 4), ' - ', round(attribute_quantile[(quant+1)], 4), "\"/>", sep=''))
+}
+qml <- append(qml, ranges)
+qml <- append(qml, "    </ranges>")
+
+#Write symbols
+qml <- append(qml, "    <symbols>")
+
+for (quant in 1:no_quantile) {
+qml <- append(qml, paste("      <symbol outputUnit=\"MM\" alpha=\"1\" type=\"line\" name=\"", (quant-1), "\">", sep=''))
+qml <- append(qml, paste("        <layer pass=\"", (quant-1), "\" class=\"SimpleLine\" locked=\"0\">", sep=''))
+qml <- append(qml, "          <prop k=\"capstyle\" v=\"square\"/>")
+qml <- append(qml, colortable[quant])
+qml <- append(qml, "          <prop k=\"customdash\" v=\"5;2\"/>")
+qml <- append(qml, "          <prop k=\"joinstyle\" v=\"bevel\"/>")
+qml <- append(qml, "          <prop k=\"offset\" v=\"0\"/>")
+qml <- append(qml, "          <prop k=\"penstyle\" v=\"solid\"/>")
+qml <- append(qml, "          <prop k=\"use_custom_dash\" v=\"0\"/>")
+qml <- append(qml, "          <prop k=\"width\" v=\"0.26\"/>")
+qml <- append(qml, "        </layer>")
+qml <- append(qml, "      </symbol>")
+}
+}
+#Write Footer
+qml <- append(qml, "    </symbols>")
+
+
+qml <- append(qml, "    <source-symbol>")
+qml <- append(qml, "      <symbol outputUnit=\"MM\" alpha=\"1\" type=\"line\" name=\"0\">")
+qml <- append(qml, "        <layer pass=\"0\" class=\"SimpleLine\" locked=\"0\">")
+qml <- append(qml, "          <prop k=\"capstyle\" v=\"square\"/>")
+qml <- append(qml, "          <prop k=\"color\" v=\"161,238,135,255\"/>")
+qml <- append(qml, "          <prop k=\"customdash\" v=\"5;2\"/>")
+qml <- append(qml, "          <prop k=\"joinstyle\" v=\"bevel\"/>")
+qml <- append(qml, "          <prop k=\"offset\" v=\"0\"/>")
+qml <- append(qml, "          <prop k=\"penstyle\" v=\"solid\"/>")
+qml <- append(qml, "          <prop k=\"use_custom_dash\" v=\"0\"/>")
+qml <- append(qml, "          <prop k=\"width\" v=\"0.26\"/>")
+qml <- append(qml, "        </layer>")
+qml <- append(qml, "      </symbol>")
+qml <- append(qml, "    </source-symbol>")
+qml <- append(qml, "    <mode name=\"quantile\"/>")
+qml <- append(qml, "    <rotation field=\"\"/>")
+qml <- append(qml, "    <sizescale field=\"\"/>")
+qml <- append(qml, "  </renderer-v2>")
+qml <- append(qml, "  <customproperties/>")
+qml <- append(qml, paste("  <displayfield>\"", attribute, "\"</displayfield>", sep=''))
+qml <- append(qml, "  <attributeactions/>")
+qml <- append(qml, "</qgis>")
+
+#Save qml-file
+qml_output <- paste(folder, paste(paste("edge_measures_", attribute, sep=''), "qml", sep='.'), sep="/")
+con_qml <- file(qml_output, open="wt")
+write.table(qml, con_qml, append = FALSE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE)
+close(con_qml)
+
+}
+
+########################################################################
+#Close R
+########################################################################
+q()
+EOF
+
+#Suppress parallel processing if not requested
+if [ $CORES -le 1 ] ; then
+	sed -i "s/mclapply/lapply/g" ${FOLDER}/network.r
+
+	sed -i "s/library(codetools)//g" ${FOLDER}/network.r
+	sed -i "s/library(multicore)//g" ${FOLDER}/network.r
+	sed -i "s/library(iterators)//g" ${FOLDER}/network.r
+	sed -i "s/library(foreach)//g" ${FOLDER}/network.r
+	sed -i "s/library(doMC)//g" ${FOLDER}/network.r
+	sed -i "s/registerDoMC()//g" ${FOLDER}/network.r
+	sed -i "s/options(cores = $CORES)//g" ${FOLDER}/network.r
+fi
+
+#Run R-igraph script
+R < ${FOLDER}/network.r --no-save --no-restore --slave
+
+#Import R-rsults to GRASS
+v.in.ogr -o --overwrite dsn="${FOLDER}/vertex_measures.vrt" output=$VERTICES
+v.in.ogr -o --overwrite dsn="${FOLDER}/edge_measures.vrt" type=line output=$EDGES
+


Property changes on: grass-addons/grass6/raster/r.connectivity.network/r.connectivity.network
___________________________________________________________________
Added: svn:executable
   + *



More information about the grass-commit mailing list