[GRASS-SVN] r65779 - in grass-addons/grass7/raster/r.green/r.green.hydro: . r.green.hydro.financial

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Jul 24 01:06:18 PDT 2015


Author: zarch
Date: 2015-07-24 01:06:18 -0700 (Fri, 24 Jul 2015)
New Revision: 65779

Added:
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/Makefile
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.html
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.py
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_compmap.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_excmap.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fcompensation.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fexcavation.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fgamma.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_flinearcost.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fnpv.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fvu.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_input.png
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_tableeco.png
Removed:
   grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/
Modified:
   grass-addons/grass7/raster/r.green/r.green.hydro/Makefile
Log:
r.green: Change r.green.hydro.economic to r.green.hydro.financial, sync with internal repo

Modified: grass-addons/grass7/raster/r.green/r.green.hydro/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/Makefile	2015-07-24 01:39:34 UTC (rev 65778)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/Makefile	2015-07-24 08:06:18 UTC (rev 65779)
@@ -5,7 +5,7 @@
 SUBDIRS = r.green.hydro.theoretical  \
           r.green.hydro.legal        \
           r.green.hydro.recommended  \
-          r.green.hydro.economic     \
+          r.green.hydro.financial     \
           r.green.hydro.closest      \
 		  r.green.hydro.optimal      \
 		  r.green.hydro.structure    \

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/Makefile	                        (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/Makefile	2015-07-24 08:06:18 UTC (rev 65779)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.financial
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/Makefile
___________________________________________________________________
Added: svn:mime-type
   + text/x-makefile
Added: svn:eol-style
   + native

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.html
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.html	                        (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.html	2015-07-24 08:06:18 UTC (rev 65779)
@@ -0,0 +1,227 @@
+<h2>DESCRIPTION</h2>
+<em>r.green.hydro.financial</em> computes the economic costs and values of the plants. It provides a cost-benefit analysis calculating realization costs and profits for each potential plant in order to know which ones are feasible.<br>
+
+<h2>NOTES</h2>
+The required input maps are the one with the segments of potentials plants (vector), the structure of these potential plants (vector), the electric grid (vector), the landuse (raster) and the slope (raster).<br><br>
+
+In the section "Input column", the column names in the table of the map with potential plants have to be reported in order to read correctly the corresponding values. <br><br>
+
+Each section of the module calculates a cost. The used formula are valid for all the currencies but the values have to be adapted. The default values related to a cost are considered in Euro. <br><br>
+
+First, we define the <b>Total cost</b> which is the sum of all the fixed costs corresponding to the construction and the implementation of the plant. It includes : 
+
+<blockquote><b>- Compensation cost :</b><br><br>
+
+That corresponds to environmental mitigation which aims at offsetting known impacts to an existing historic or natural resource such as a stream, wetland, endangered species, archeological site or historic structure. Environmental mitigation is typically a part of an environmental crediting system established by governing bodies which involves allocating debits and credits. Debits occur in situations where a natural resource has been destroyed or severely impaired. When an entity such as a business or individual has a "debit" they are required to purchase a "credit" to mitigate (offset) the damages they did. In our case, the compensation cost can be calculated thanks to this formula :<br><br>
+<center>
+<img src="r_green_hydro_financial_fcompensation.png" alt="formula_compensation"><br></center><br>
+<blockquote>where L<SUB>v</SUB> is a raster map with the land use value [currency/ha];<br>
+T<SUB>r</SUB> is a raster map with the tributes [currency/ha];<br>
+V<SUB>u</SUB> is a raster map with the value upper part of the soil [currency/ha];<br>
+r is a scalar with the interest rate (default value: 0.03);<br>
+n is the life of hydropower plant [year] (default value: 30);<br>
+γ is a scalar (default value: 5/4);<br>
+w<SUB>e</SUB> is a scalar with the average width excavation [m] (default value: 2);<br>
+res is a scalar with the raster resolution.<br><br>
+
+The V<SUB>u</SUB> raster map is computed with the formula:<br><br>
+<center>
+<img src="r_green_hydro_financial_fvu.png" alt="Vu"><br></center><br>
+<blockquote> where S<SUB>v</SUB> is a raster map with the stumpage value [currency/ha];<br>
+Rot is a raster map with the value rotation period per land use type [year];<br>
+Y is a raster map with the average year [year];</blockquote><br>
+
+The user can add directly the maps L<SUB>v</SUB>, T<SUB>r</SUB>, S<SUB>v</SUB>, Rot and Y as inputs. Otherwise, the maps can be computed using the land use raster map and reclassifying the values with the module r.reclass. The program creates the reclassified maps if the user provides the input text files for each category (the input data is the path of the text file). Here is an example of a text file to create the landvalue raster map, the costs are in currency/ha :<br>
+<pre>1 = 0 rocks, macerated, glaciers
+2 = 0 urbanized areas, infrastructure
+3 = 0 shores
+4 = 0 waters
+5 = 200 gardens
+6 = 4000 mining areas
+7 = 2000 agricultural areas
+8 = 1500 meadows
+9 = 1000 areas with predominantly pastoral value
+10 = 3000  forestry land</pre></blockquote>
+
+Once the calculation is made, a new column with the compensation cost is added in the table of the input map with potential plants. A raster map with the compensation cost can also be computed, as well as a raster map with the value upper part of the soil (See Optional section).<br><br><br>
+
+<b>- Excavation cost :</b><br><br>
+
+This cost concerns the works of digging to implement channels. It is calculated as followed :<br><br>
+<center><img src="r_green_hydro_financial_fexcavation.png" alt="Excavation"><br></center><br>
+<blockquote>where S is a raster map with the slope in [%];<br>
+Ψ is a raster map with values of minimum excavation costs [currency/mc]<br>
+Φ is a raster map with values of maximum excavation costs [currency/mc]<br>
+w is the width of the excavation [m] (default value: 2);<br>
+d is the depth of the excavation [m] (default value: 2);<br>
+l is the length of the excavation [m] which depends on the channels lengths;<br><br></blockquote>
+
+If the user hasn't the raster maps Ψ and Φ, the latter can be computed from the land use raster map if the user provides a text file with the reclassification values (from land use value to excavation cost (min or max)). It is the same principle as the one explained above for landvalue, tributes, stumpage, rotation period per land use type and average age.<br><br>
+
+The user can choose to put a slope limit (S<SUB>lim</SUB>) [%] above which the cost will be equal to the maximum cost.<br><br>
+
+Then, a new column with the excavation cost is added in the table of the input map with potential plants. A raster map with excavation cost can also be computed (See Optional section).<br><br><br>
+
+<b>- Electro-mechanical cost :</b><br><br>
+
+It is the cost of the electro-mechanical equipment which includes the turbine, the alternator and the regulator. It is a high percentage of a small hydropower plant budget (around 30% and 40% of the total sum). <br><br>
+
+We use the Aggidis et al. formula which calculates this cost thanks to the values of power and head :<br><br>
+<center>CEM = γ * power<SUP>α</SUP> * head<SUP>β</SUP> + c</center>
+<blockquote> where power is the power of the plant [kW];<br>
+head is the head of the plant [m];<br>
+α is a power coefficient (default value: 0.56);<br>
+β is a head coefficient (default value: -0.112);<br>
+γ is a coefficient (default value: 15600);<br>
+c is a constant (default value: 0);<br></blockquote>
+
+A new column with the electro-mechanical costs is added in the table of the input map with potential plants.<br><br><br>
+
+<b>- Supply and installation cost for pipeline and electroline :</b><br><br>
+
+This is the sum of the supply and installation costs for the derivation channel, the penstock (boths compose the pipeline) and the electroline which links the transformer near the turbine to the existing grid. The formula is the following :<br><br>
+<center><img src="r_green_hydro_financial_flinearcost.png" alt="Lines cost"><br></center>
+<blockquote>where α is the supply and installation costs [currency/m] (default value for pipeline: 310 Euro/m, and for the electroline 250 Euro/m);<br>
+l is the length of the line [m], the pipeline length is found in the structure map and the electroline length is computed thanks to the grid vector map;<br><br><br></blockquote>
+
+
+<b>- Power station cost</b><br><br>
+
+It concerns the construction cost of the building composing the power station. It is considered as a percentage of the electro-mechanical cost :<br><br>
+<center>CST = α * CEM</center>
+<blockquote>where α is as default 0.52;</blockquote><br><br><br>
+
+
+<b>- Inlet cost</b><br><br>
+
+It concerns the construction cost of the water intake structure. It is considered as a percentage of the electro-mechanical cost :<br><br>
+<center>CIN = α *EM</center>
+<blockquote>where α is as default 0.38;</blockquote></blockquote><br><br>
+
+
+All these costs define the <b>Total cost</b> :<br><br>
+<center>Total cost = (CCM+CL+CEX+CEM+CST+CIN+CGRD)(1+α+β)</center>
+<blockquote>where CCM is the Compensation cost [currency];<br>
+CL is the Supply and installation cost for pipeline and electroline [currency];<br>
+CEX is the Excavation cost [currency];<br>
+CEM is the Electro-mechanical cost [currency];<br>
+CST is the Power station cost [currency];<br>
+CIN is the Inlet cost [currency];<br>
+CGRD is the Grid connection cost, (default is 50000 Euro), that includes the easement indemnity;<br>
+α is the factor to consider general expenses, (default is 0.15);<br>
+β is the factor to consider hindrances expenses, (default is 0.10);<br><br></blockquote>
+
+α and β are factors which offset the underestimation of this total cost. Indeed, some other costs have to be taken into consideration but it's hard to make a general assessment because they are specific for each plant. Moreover, for each plant realization there are unexpected events (hindrances) which make the implementation more complex and expensive.<br><br><br>
+
+Then, the module calculates the <b>maintenance cost per year</b>, according to this formula :<br><br>
+<center>maintenance = α * power<SUP>1+β</SUP> + c </center>
+<blockquote>where power is the power of the plant [kW]; <br>
+α is a power coefficient (default value: 3871.2); <br>
+β is a power coefficient (default value: 0.45); <br>
+c is a constant, (default value: 0); <br><br><br></blockquote>
+
+The <b>yearly revenue</b> corresponds to the money gained selling all the electricity the plant produces in a year. It is simply calculated as the product of the power produced in a year by the price of the electricity (including some coefficients) :<br><br>
+<center>revenue = η * power * price * yh * α + c </center>
+<blockquote> where η is the efficiency of the electro-mechanical components (default value: 0.81); <br>
+power is the installed power of the plant [kW]; <br>
+price is the energy price [currency/kW], (default value: 0.09 Euro/kW); <br>
+yh is the yearly operation hours of the plant [hours], (default value: 6500); <br>
+α is the coefficient to transform installed power to mean power, (default value: 
+0.5); <br>
+c is a constant, (default value: 0);</blockquote><br><br>
+
+Finally, all these values allow to calculate the <b>Net Present Value (NPV)</b>. It is the sum of the present values of incoming and outgoing cash flows over a period of time. It allow to know if there are profits so if the plant is feasible. It corresponds to :<br><br>
+<center><img src="r_green_hydro_financial_fnpv.png" alt="NPV"><br></center>
+<blockquote>where revenue is the yearly revenue value [currency/year]; <br>
+maintenance is the yearly maintenance value [currency/year]; <br>
+Total cost is the total cost of the plant [currency] <br>
+n the number of years of the plant [year] (default value: 30) <br>
+γ is a coefficient which assesses the cost of interest and amortization. It is defined as:<br><br>
+<center><img src="r_green_hydro_financial_fgamma.png" alt="Gamma"><br></center><br>
+<blockquote>where r is the interest rate (default value: 0.03).<br><br></blockquote>
+</blockquote>
+
+More concretely, the program computes the following results :<br>
+<blockquote>- the input map with the structure of the plants has an updated table with the different costs of construction and implementation and their sum (tot_cost).<br><br>
+
+- the created output map shows the structure of the potential plants with a re-organized table. The latter doesn't make the difference between derivation channel and penstock. Each line gives the intake_id, plant_id, side (structures are computed on both sides of the river), power (kW), gross_head (m), discharge (m<SUP>3</SUP>/s), tot_cost (total cost for construction and implementation), yearly maintenance cost, yearly revenue,net present value (NPV) and max_NPV. The structure of potential plants is given for each side of the river, max_NPV is 'yes' for the side with the highest NPV and 'no' for the other side.<br><br>
+
+- the input map with the segments of the plants has an updated table with the total cost, yearly maintenance cost, yearly revenue and the net present value. The parameter "segment_basename" (in Input column) allows to add a prefix at the column names in order to show the results for different cases in the same table without overwriting the columns.<br><br>
+
+- In the Optional section, there is the possibility to create three raster maps showing the compensation, excavation and upper part of the soil values.</blockquote><br>
+
+<h2>EXAMPLE</h2>
+
+This example is based on the case-study of Gesso and Vermenagna valleys in the Natural Park of the Maritime Alps, Piedmont, Italy.<br><br>
+
+Here in black is the input vector file techplants with the structure of the potential plants and the technical value of power (including head losses and efficiencies, computed by r.green.hydro.technical). The vector map with the segments of river potentialplants is also visibile in blue on this picture, as well as the vector map with the intakes and restitutions of the plants (in red).
+
+<center>
+<img src="r_green_hydro_financial_input.png" alt="input"><br>
+Structure of the potential plants
+</center><br><br>
+
+The following command updates the table of structplants and segplants adding the costs :<br>
+<div class="code"><pre>
+r.green.hydro.financial 
+plant=potentialplants 
+struct=techplants
+plant_head_column=net_head 
+landuse=landuse 
+rules_landvalue=/pathtothefile/landvalue.rules 
+rules_tributes=pathtothefile/tributes.rules 
+rules_stumpage=/pathtothefile/stumpage.rules 
+rules_rotation=/pathtothefile/rotation.rules 
+rules_age=/pathtothefile/age.rules 
+slope=slope 
+rules_min_exc=/pathtothefile/excmin.rules 
+rules_max_exc=/pathtothefile/excmax.rules 
+electro=grid 
+output_struct=ecoplants 
+compensation=comp 
+excavation=exc 
+upper=upper</pre></div><br>
+
+It also creates four new raster maps (ecoplants, comp, exc and upper) :<br>
+<blockquote>- ecoplants which shows the structure of the potential plants. The table contains these four columns (total cost, maintenance cost, revenue and NPV) :<br><br>
+<center>
+<img src="r_green_hydro_financial_tableeco.png" alt="table_eco"><br>
+Table of the ouput raster map ecoplants
+</center><br>
+
+The same columns are added for the input map with the segments (potentialplants). In the table of the input map with the structure (techplants), the different costs which compose the total cost are added in columns (but not the four previous values).<br><br>
+
+- comp which shows the compensation values (in currency) for each land use :<br>
+<center>
+<img src="r_green_hydro_financial_compmap.png" alt="Compensation_map"><br>
+Output raster map with compensation values
+</center><br><br>
+- upper which shows the values of the upper part of the soil (in currency) for each land use.<br><br>
+- exc which shows the excavation value (in currency) for each land use.<br><br>
+
+These two last maps look like comp, but with their corresponding values.<br><br></blockquote>
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="r.green.hydro.discharge.html">r.green.hydro.discharge</a><br>
+<a href="r.green.hydro.delplants.html">r.green.hydro.delplants</a><br>
+<a href="r.green.hydro.theoretical.html">r.green.hydro.theoretical</a><br>
+<a href="r.green.hydro.optimal.html">r.green.hydro.optimal</a><br>
+<a href="r.green.hydro.recommended.html">r.green.hydro.recommended</a><br>
+<a href="r.green.hydro.structure.html">r.green.hydro.structure</a><br>
+<a href="r.green.hydro.technical.html">r.green.hydro.technical</a><br>
+</em>
+
+<h2>REFERENCE</h2>
+<i>The costs of small-scale hydro power production: Impact on the development
+of existing potential</i>, p.2635, Aggidis et Al, 2010<br><br>
+Text on expropriation for public utility D.P.R. n. 327/2001, updated in 2013<br><br>
+The excavation costs are found in the price list of the Piedmont Region for public works<br><br>
+The costs of intake and civil works are based on analysis of technical documents for mini-hydro projects in Italy<br><br>
+The calculation of the yearly maintenance cost is based on the report AEEG for the Polytechnic University of Milan<br><br>
+The definition of environmental mitigation comes from Wikipedia
+
+
+<h2>AUTHORS</h2>
+Sandro Sacchelli (University of Florence, Italy), Pietro Zambelli (Eurac Research, Bolzano, Italy), Manual written by Julie Gros.<br>
+Last changed: $Date : 2015-07-07 16:13 GMT+1$


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.html
___________________________________________________________________
Added: svn:mime-type
   + text/html
Added: svn:keywords
   + Author Date Id
Added: svn:eol-style
   + native

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.py	                        (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.py	2015-07-24 08:06:18 UTC (rev 65779)
@@ -0,0 +1,1043 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE:      r.green.hydro.financial
+# AUTHOR(S):   Sandro Sacchelli (BASH-concept), Pietro Zambelli (Python)
+# PURPOSE:     Assess the hydro plant costs
+# COPYRIGHT:   (C) 2014 by the GRASS Development Team
+#
+#              This program is free software under the GNU General Public
+#              License (>=v2). Read the file COPYING that comes with GRASS
+#              for details.
+#
+#############################################################################
+#
+#%Module
+#% description: Assess the financial costs and values
+#% overwrite: yes
+#%End
+#%option G_OPT_V_INPUT
+#%  key: plant
+#%  label: Name of the input vector map with the segments of the plants
+#%  required: yes
+#%end
+#%option G_OPT_V_INPUT
+#%  key: struct
+#%  label: Name of the input vector map with the structure of the plants
+#%  required: yes
+#%end
+
+#############################################################################
+# DEFINE COLUMNS OF V INPUT
+
+
+# TODO: change the order as disccussed with Giulia and Francesco
+# the order should be: {vector name}_{property}_{extra}
+# with:
+# * vector name the string use to identify the vector map
+# * property: layer, column
+# * other string to clarify the key
+#
+# therefore:
+# * plant_id_column => plant_column_id
+# * plant_power_column => plant_column_power
+# etc.
+
+#%option G_OPT_V_FIELD
+#%  key: segment_layer
+#%  label: Name of the vector map layer of the segments
+#%  required: no
+#%  answer: 1
+#%  guisection: Input columns
+#%end
+#%option G_OPT_V_FIELD
+#%  key: plant_layer
+#%  label: Name of the vector map layer of the structure of the plants
+#%  required: no
+#%  answer: 1
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_id_column
+#%  type: string
+#%  description: Column name with the plant id
+#%  required: no
+#%  answer: plant_id
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_power_column
+#%  type: string
+#%  description: Column name with power value
+#%  required: no
+#%  answer: power
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_head_column
+#%  type: string
+#%  description: Column name with head value
+#%  required: no
+#%  answer: gross_head
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_side_column
+#%  type: string
+#%  description: Column name with the strings that define the side of the plant
+#%  required: no
+#%  answer: side
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_kind_column
+#%  type: string
+#%  description: Column name with the strings that define if it's a derivation channel or a penstock
+#%  required: no
+#%  answer: kind
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_kind_intake
+#%  type: string
+#%  description: Value contained in the column 'kind' which corresponds to the derivation channel
+#%  required: no
+#%  answer: conduct
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: plant_kind_turbine
+#%  type: string
+#%  description: Value contained in the column 'kind' which corresponds to the penstock
+#%  required: no
+#%  answer: penstock
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: segment_column_id
+#%  description: Optional map with the segments : column name with the plant id
+#%  required: no
+#%  answer: plant_id
+#%  guisection: Input columns
+#%end
+#%option
+#%  key: segment_basename
+#%  type: string
+#%  description:Optional map with the segments : basename of the columns that will be added to the segment vector map
+#%  required: no
+#%  answer: case1
+#%  guisection: Input columns
+#%end
+
+
+#############################################################################
+# DEFINE COMPENSATION COSTS
+# provide raster
+
+#%option
+#%  key: interest_rate
+#%  type: double
+#%  description: Interest rate value
+#%  required: no
+#%  answer: 0.03
+#%  guisection: Compensation
+#%end
+#%option
+#%  key: gamma_comp
+#%  type: double
+#%  description: Coefficient
+#%  required: no
+#%  answer: 1.25
+#%  guisection: Compensation
+#%end
+#%option
+#%  key: life
+#%  type: double
+#%  description: Life of the hydropower plant [year]
+#%  required: no
+#%  answer: 30
+#%  guisection: Compensation
+#%end
+
+
+#%option G_OPT_R_INPUT
+#%  key: landvalue
+#%  label: Name of the raster map with the land value [currency/ha]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#%  key: tributes
+#%  label: Name of the raster map with the tributes [currency/ha]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#%  key: stumpage
+#%  label: Name of the raster map with the stumpage value [currency/ha]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#%  key: rotation
+#%  label: Name of the raster map with the rotation period per landuse type [year]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#%  key: age
+#%  label: Name of the raster map with the average age [year]
+#%  required: no
+#%  guisection: Compensation
+#%end
+
+# or rule to transform landuse categories in raster maps
+#%option G_OPT_R_INPUT
+#%  key: landuse
+#%  label: Name of the raster map with the landuse categories
+#%  required: no
+#%  guisection: Compensation
+#%end
+
+# RULES
+#%option
+#%  key: rules_landvalue
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate a value for each landuse [currency/ha]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option
+#%  key: rules_tributes
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate a tribute rate for each landuse [currency/ha]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option
+#%  key: rules_stumpage
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate a stumpage value for each landuse [currency/ha]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option
+#%  key: rules_rotation
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate a rotation value for each landuse [year]
+#%  required: no
+#%  guisection: Compensation
+#%end
+#%option
+#%  key: rules_age
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate an age for each landuse [year]
+#%  required: no
+#%  guisection: Compensation
+#%end
+
+#############################################################################
+# DEFINE EXCAVATION COSTS
+
+#%option
+#%  key: width
+#%  type: double
+#%  description: Width of the excavation works [m]
+#%  required: no
+#%  answer: 2.
+#%  guisection: Excavation
+#%end
+#%option
+#%  key: depth
+#%  type: double
+#%  description:Depth of the excavation works [m]
+#%  required: no
+#%  answer: 2.
+#%  guisection: Excavation
+#%end
+#%option
+#%  key: slope_limit
+#%  type: double
+#%  description: Slope limit, above this limit the cost will be equal to the maximum [degree]
+#%  required: no
+#%  answer: 50.
+#%  guisection: Excavation
+#%end
+
+# provide raster maps
+
+#%option G_OPT_R_INPUT
+#%  key: min_exc
+#%  label: Minimum excavation costs [currency/mc]
+#%  required: no
+#%  guisection: Excavation
+#%end
+#%option G_OPT_R_INPUT
+#%  key: max_exc
+#%  label: Maximum excavation costs [currency/mc]
+#%  required: no
+#%  guisection: Excavation
+#%end
+#%option G_OPT_R_INPUT
+#%  key: slope
+#%  label: Slope raster map
+#%  required: yes
+#%  guisection: Excavation
+#%end
+
+
+
+# RULES
+#%option
+#%  key: rules_min_exc
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate a minimum excavation cost for each landuse [currency/mc].
+#%  required: no
+#%  guisection: Excavation
+#%end
+#%option
+#%  key: rules_max_exc
+#%  type: string
+#%  description: Rule file for the reclassification of the landuse to associate a maximum excavation cost for each landuse [currency/mc].
+#%  required: no
+#%  guisection: Excavation
+#%end
+
+#############################################################################
+# DEFINE ELECTROMECHANICAL COSTS
+#%option
+#%  key: alpha_em
+#%  type: double
+#%  description: Electro-mechanical costs alpha parameter, default values taken from Aggidis et al. 2010
+#%  required: no
+#%  answer: 0.56
+#%  guisection: Electro-mechanical
+#%end
+#%option
+#%  key: beta_em
+#%  type: double
+#%  description: Electro-mechanical costs beta parameter, default values taken from Aggidis et al. 2010
+#%  required: no
+#%  answer: -0.112
+#%  guisection: Electro-mechanical
+#%end
+#%option
+#%  key: gamma_em
+#%  type: double
+#%  description: Electro-mechanical costs gamma parameter, default values taken from Aggidis et al. 2010
+#%  required: no
+#%  answer: 15600.0
+#%  guisection: Electro-mechanical
+#%end
+#%option
+#%  key: const_em
+#%  type: double
+#%  description: Electro-mechanical costs constant value, default values taken from Aggidis et al. 2010
+#%  required: no
+#%  answer: 0.
+#%  guisection: Electro-mechanical
+#%end
+
+#############################################################################
+# DEFINE SUPPLY AND INSTALLATION COSTS
+#%option
+#%  key: lc_pipe
+#%  type: double
+#%  description: Supply and installation linear cost for the pipeline [currency/m]
+#%  answer: 310.
+#%  guisection: Supply & Installation
+#%end
+#%option
+#%  key: lc_electro
+#%  type: double
+#%  description: Supply and installation linear cost for the electroline [currency/m]
+#%  answer: 250.
+#%  guisection: Supply & Installation
+#%end
+#%option G_OPT_V_INPUT
+#%  key: electro
+#%  label: Name of the vector map with the electric grid
+#%  required: yes
+#%end
+#%option
+#%  key: electro_layer
+#%  description: Name of the vector map layer of the grid
+#%  required: no
+#%  answer: 1
+#%end
+
+#############################################################################
+# DEFINE POWER STATION COSTS
+#%option
+#%  key: alpha_station
+#%  type: double
+#%  description: Power station costs are assessed as a fraction of the Electro-mechanical costs
+#%  answer: 0.52
+#%  guisection: Power station
+#%end
+
+#############################################################################
+# DEFINE INLET COSTS
+#%option
+#%  key: alpha_inlet
+#%  type: double
+#%  description: Inlet costs are assessed as a fraction of the Electro-mechanical costs
+#%  answer: 0.38
+#%  guisection: Inlet
+#%end
+
+#############################################################################
+# DEFINE OTHER COSTS
+#%option
+#%  key: grid
+#%  type: double
+#%  description:  Cost for grid connection
+#%  answer: 50000
+#%  guisection: Other
+#%end
+#%option
+#%  key: general
+#%  type: double
+#%  description:  Factor for general expenses
+#%  answer: 0.15
+#%  guisection: Other
+#%end
+#%option
+#%  key: hindrances
+#%  type: double
+#%  description:  Factor for hindrances expenses
+#%  answer: 0.1
+#%  guisection: Other
+#%end
+# DEFINE MAINTENANCE COSTS
+#%option
+#%  key: cost_maintenance_per_kw
+#%  type: double
+#%  description: Maintenace costs per kW
+#%  answer: 7000.
+#%  guisection: Maintenance
+#%end
+#%option
+#%  key: alpha_maintenance
+#%  type: double
+#%  description: Alpha coefficient to assess the maintenance costs
+#%  answer: 0.05
+#%  guisection: Maintenance
+#%end
+#%option
+#%  key: beta_maintenance
+#%  type: double
+#%  description: Beta coefficient to assess the maintenance costs
+#%  answer: 0.45
+#%  guisection: Maintenance
+#%end
+#%option
+#%  key: const_maintenance
+#%  type: double
+#%  description: Constant to assess the maintenance costs
+#%  answer: 0.
+#%  guisection: Maintenance
+#%end
+
+# DEFINE REVENUES PARAMETERS
+#%option
+#%  key: energy_price
+#%  type: double
+#%  description: Energy price per kW [currency/kW]
+#%  answer: 0.1
+#%  guisection: Revenues
+#%end
+#%option
+#%  key: eta
+#%  type: double
+#%  description: Efficiency of electro-mechanical components
+#%  answer: 0.81
+#%  guisection: Revenues
+#%end
+#%option
+#%  key: operative_hours
+#%  type: double
+#%  description: Number of operative hours per year [hours/year]
+#%  answer: 6500.
+#%  guisection: Revenues
+#%end
+#%option
+#%  key: alpha_revenue
+#%  type: double
+#%  description: Coefficient to transform installed power to mean power
+#%  answer: 1.0
+#%  guisection: Revenues
+#%end
+#%option
+#%  key: const_revenue
+#%  type: double
+#%  description: Constant to assess the revenues
+#%  answer: 0.
+#%  guisection: Revenues
+#%end
+
+#############################################################################
+# DEFINE OUTPUTS
+#%option G_OPT_V_OUTPUT
+#%  key: output_struct
+#%  label: Name of the output vector map : plants' structure including the main costs in the table
+#%  required: yes
+#%end
+
+## COSTS
+#%option G_OPT_R_OUTPUT
+#%  key: compensation
+#%  label: Raster map with the compensation values
+#%  required: no
+#%end
+#%option G_OPT_R_OUTPUT
+#%  key: excavation
+#%  label: Raster map with the excavation costs
+#%  required: no
+#%end
+
+## VALUES
+#%option G_OPT_R_OUTPUT
+#%  key: upper
+#%  label: Raster map with the value upper part of the soil 
+#%  required: no
+#%end
+
+#############################################################################
+from __future__ import print_function
+
+import os
+import atexit
+import random
+
+import numpy as np
+
+from grass.exceptions import ParameterError
+from grass.script.core import parser, overwrite, warning
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+
+from grass.pygrass.utils import set_path
+from grass.pygrass.raster import RasterRow
+from grass.pygrass.vector import VectorTopo, sql
+from grass.pygrass.vector.basic import Cats
+
+from grass.pygrass.messages import get_msgr
+
+try:
+    import numexpr as ne
+except ImportError:
+    warning('You should install numexpr to use this module: '
+                  'pip install numexpr')
+
+# set python path to the shared r.green libraries
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+
+from libgreen.utils import cleanup
+from libgreen.checkparameter import check_required_columns, exception2error
+from libhydro.plant import read_plants, write_structures
+
+
+def rname(base):
+    return '%s_%03d' % (base, random.randint(0, 1000))
+
+
+def check_raster_or_landuse(opts, params):
+    """Return a list of raster name, i f necessary the rasters are generated
+    using the file with the rules to convert a landuse map to the raster
+    that it is needed.
+    """
+    msgr = get_msgr()
+    rasters = []
+    for par in params:
+        if opts[par]:
+            rasters.append(opts[par])
+        elif opts['rules_%s' % par] and opts['landuse']:
+            output = rname(par)
+            msg = 'Creating: {out} using the rules: {rul}'
+            msgr.verbose(msg.format(out=output, rul='rules_%s' % par))
+            r.reclass(input=opts['landuse'], output=output,
+                      rules=opts['rules_%s' % par])
+            rasters.append(output)
+        else:
+            msg = '{par} or rule_{par} are required'
+            raise ParameterError(msg.format(par=par))
+    return rasters
+
+
+def upper_value(upper, stu, lan, rot, age, irate, overwrite=False):
+    """Compute the upper value of a land.
+    """
+    expr = ("{upper} = ({stu} + {lan}) / ((1 + {irate})^({rot} - {age}) )"
+            " - {lan}")
+    r.mapcalc(expr.format(upper=upper, stu=stu, lan=lan, irate=irate,
+                          rot=rot, age=age), overwrite=overwrite)
+
+
+def compensation_cost(comp, lan, tri, upper,
+                      irate, gamma, life, width, overwrite=False):
+    """Compute the compensation raster map costs"""
+    expr = ("{comp} = (({lan} + {tri} * "
+            "(1 + {irate})^{life}/({irate} * (1 + {irate}))) "
+            "* {gamma} + {upper}) * {width} * nsres() / 10000")
+    r.mapcalc(expr.format(comp=comp, lan=lan, tri=tri, upper=upper,
+                          irate=irate, gamma=gamma, life=life, width=width),
+              overwrite=overwrite)
+
+
+def excavation_cost(exc, excmin, excmax, slope,
+                    slim, width, depth, overwrite=False):
+    """Compute the excavation cost"""
+    expr = ("{exc} = if({slope} < {slim}, "
+            "({excmin} + ({excmax} - {excmin}) / {slim} * {slope})"
+            "* {width} * {depth} * nsres(),"
+            "{excmax} * {width} * {depth} * nsres())")
+    r.mapcalc(expr.format(exc=exc, slope=slope, excmin=excmin, excmax=excmax,
+                          slim=slim, width=width, depth=depth),
+              overwrite=overwrite)
+
+def vmapcalc2(vmap, vlayer, cname, ctype, expr, overwrite=False):
+    v.db_addcolumn(map=vmap, layer=vlayer, columns=[(cname, ctype)],
+                   overwrite=overwrite)
+    v.db_update(map=vmap, layer=vlayer, column=cname, query_column=expr,
+                overwrite=overwrite)
+
+
+def get_cnames(expr,
+               _names_cache=ne.utils.CacheDict(256),
+               _numexpr_cache=ne.utils.CacheDict(256), **kwargs):
+    if not isinstance(expr, (str, unicode)):
+        raise ValueError("must specify expression as a string")
+    # Get the names for this expression
+    context = ne.necompiler.getContext(kwargs, frame_depth=1)
+    expr_key = (expr, tuple(sorted(context.items())))
+    if expr_key not in _names_cache:
+        _names_cache[expr_key] = ne.necompiler.getExprNames(expr.strip(),
+                                                            context)
+    names, ex_uses_vml = _names_cache[expr_key]
+    return names
+
+
+def vcolcalc(vname, vlayer, ctype, expr, condition=lambda x: x is None,
+             notfinitesubstitute=None, **kwargs):
+    equal = expr.index('=')
+    if equal < 0:
+        raise
+    cname = expr[:equal].strip()
+    expr = expr[(equal + 1):].strip()
+    cnames = get_cnames(expr)
+    vname, vmapset = vname.split('@') if '@' in vname else (vname, '')
+    with VectorTopo(vname, mapset=vmapset, layer=vlayer, mode='r') as vect:
+        if vect.table is None:
+            msg = 'Vector: {vname} is without table.'
+            raise TypeError(msg.format(vname=vname))
+
+        cols = vect.table.columns
+        # check if the column in the expressions exist or not
+        for col in cnames:
+            if col not in cols:
+                msg = 'Vector: {vname} has not column: {col} in layer: {layer}'
+                raise  TypeError(msg.format(vname=vname, col=col,
+                                            layer=vlayer))
+        if cname not in cols:
+            vect.table.columns.add(cname, ctype)
+        # extract value from the attribute table
+        if cnames[0] != vect.table.key:
+            cnames.insert(0, vect.table.key)
+        # TODO: find a more elegant way to do it
+        sql = vect.table.filters.select(', '.join(cnames)).where(' is not NULL AND '.join(cnames)).get_sql()[:-1] + ' is not NULL;'
+        data = np.array(list(vect.table.execute(sql)))
+        # create a dictionary with local variables
+        lvars = {col: array for col, array in zip(cnames, data.T)}
+        # compute the result with numexpr
+        res = ne.evaluate(expr, local_dict=lvars, **kwargs)
+        if notfinitesubstitute is not None:
+            res[~np.isfinite(res)] = notfinitesubstitute
+        # save the result to the vector point
+        cur = vect.table.conn.cursor()
+        upstr = 'UPDATE {tname} SET {cname}=? WHERE {key} == ?;'
+        cur.executemany(upstr.format(tname=vect.table.name, cname=cname,
+                                     key=vect.table.key),
+                        zip(res, lvars[vect.table.key]))
+        # save changes
+        vect.table.conn.commit()
+
+
+def electromechanical_cost(vname, power, head,
+                           gamma=15600., alpha=0.56, beta=0.112, const=0.,
+                           vlayer=1, cname='em_cost',
+                           ctype='double precision',
+                           overwrite=False):
+    expr = "{cname} = {gamma} * {power}**{alpha} * {head}**{beta} + {const}"
+    vcolcalc(vname, vlayer, ctype, notfinitesubstitute=0.,
+             expr=expr.format(cname=cname, gamma=gamma, power=power,
+                              alpha=alpha, head=head, beta=beta, const=const))
+
+
+def col_exist(vname, cname, ctype='double precision', vlayer=1, create=False):
+    vname, vmapset = vname.split('@') if '@' in vname else (vname, '')
+    with VectorTopo(vname, mapset=vmapset, layer=vlayer, mode='r') as vect:
+        if vect.table is None:
+            msg = 'Vector: {vname} is without table.'
+            raise TypeError(msg.format(vname=vname))
+        res = cname in vect.table.columns
+        if res is False:
+            vect.table.columns.add(cname, ctype)
+        return res
+
+
+def linear_cost(vname, cname='lin_cost', alpha=310., length='length', vlayer=1,
+                ctype='double precision',  overwrite=False):
+    # check if length it is alread in the db
+    if not col_exist(vname, 'length', create=True):
+        v.to_db(map=vname, type='line', layer=vlayer, option='length',
+                columns='length')
+    expr = "{cname} = {alpha} * {length}"
+    vcolcalc(vname, vlayer, ctype, notfinitesubstitute=0.,
+             expr=expr.format(cname=cname, alpha=alpha, length=length))
+
+
+def get_electro_length(opts):
+    # open vector plant
+    pname = opts['struct']
+    pname, vmapset = pname.split('@') if '@' in pname else (pname, '')
+    with VectorTopo(pname, mapset=vmapset, layer=int(opts['plant_layer']),
+                    mode='r') as vect:
+        kcol = opts['plant_kind_column']
+        ktype = opts['plant_kind_turbine']
+        # check if electro_length it is alredy in the table
+        if 'electro_length' not in vect.table.columns:
+            vect.table.columns.add('electro_length', 'double precision')
+        # open vector map with the existing electroline
+        ename = opts['electro']
+        ename, emapset = ename.split('@') if '@' in ename else (ename, '')
+        with VectorTopo(ename, mapset=emapset, layer=int(opts['electro_layer']),
+                        mode='r') as electro:
+            for line in vect:
+                if line.attrs[kcol] == ktype:
+                    # the turbine is the last point of the penstock
+                    turbine = line[-1]
+                    # find the closest electro line
+                    eline = electro.find['by_point'].geo(turbine, maxdist=1e6)
+                    dist = eline.distance(turbine)
+                    line.attrs['electro_length'] = dist.dist
+                else:
+                    line.attrs['electro_length'] = 0.
+            vect.table.conn.commit()
+
+
+def get_gamma_NPV(r=0.03, y=30):
+    """gamma it is a coefficient define as:
+
+        $\gamma = 1 - \frac{1 - (1+r)^y}{r(1+r)^y}$
+
+    with $r$ as the interest rate (default value: 0.03) and
+    $y$ as the number of years of the plant [years] (default value: 30);
+    """
+    return 1 - (1 - (1 + r)**y) / (r * (1 + r)**y)
+
+
+def group_by(vinput, voutput, isolate=None, aggregate=None,
+             function='sum', vtype='lines',
+             where='',  group_by=None, linput=1, loutput=1):
+    vname, vmapset = vinput.split('@') if '@' in vinput else (vinput, '')
+    with VectorTopo(vname, mapset=vmapset, mode='r') as vin:
+        columns = ['cat', ]
+        vincols = vin.table.columns
+        types = ['PRIMARY KEY', ]
+        siso = ''
+        if isolate is not None:
+            columns += list(isolate)
+            types += [vincols[c] for c in isolate]
+            siso = ', '.join(isolate)
+        sagg = ''
+        if aggregate is not None:
+            ct = '{func}({col}) as {col}'
+            sagg = ', '.join([ct.format(func=function, col=col)
+                              for col in aggregate])
+            columns += list(aggregate)
+            types += [vincols[c] for c in aggregate]
+
+        scols = "%s, %s" % (siso, sagg) if siso and sagg else siso + sagg
+        base = "SELECT {cols} FROM {tin}"
+        bwhere = " WHERE %s" % where if where else ''
+        bgroup = " GROUP BY %s" % ', '.join(group_by) if group_by else ''
+        grp = (base + bwhere + bgroup + ';').format(cols=scols,
+                                                    tin=vin.table.name,
+                                                    tout=voutput)
+        bqry = "SELECT {cat} FROM {tin} WHERE {cond};"
+        qry = bqry.format(cat=vin.table.key, tin=vin.table.name,
+                          cond=' AND '.join(['%s=?' % g for g in group_by]))
+        selcols = columns[1:]
+        gindexs = [selcols.index(col) for col in group_by]
+        insrt = sql.INSERT.format(tname=voutput,
+                                  values=','.join(['?', ] * len(columns)))
+        with VectorTopo(voutput, mode='w', tab_name=voutput, layer=loutput,
+                        tab_cols=list(zip(columns, types)), link_key='cat',
+                        overwrite=True) as vout:
+            cur = vout.table.conn.cursor()
+            ncat = 1
+            #import ipdb; ipdb.set_trace()
+            # TODO: why do I need to use list(cur) and I can not iterate directly
+            for row in list(cur.execute(grp)):
+                # add the new line to table
+                print(row)
+                cur.execute(insrt, tuple([ncat, ] + list(row)))
+                for cat in cur.execute(qry, tuple([row[i] for i in gindexs])):
+                    for line in vin.cat(cat[0], vtype):
+                        # set the new category
+                        category = Cats(line.c_cats)
+                        category.reset()
+                        category.set(ncat, loutput)
+                        # write geometry
+                        vout.write(line, set_cats=False)
+                ncat += 1
+            vout.table.conn.commit()
+
+
+def max_NPV(l0, l1):
+    return l0 if l0.attrs['NPV'] > l1.attrs['NPV'] else l1
+
+
+def economic2segment(economic, segment, basename='eco_',
+                     eco_layer=1, seg_layer=1,
+                     eco_pid='plant_id', seg_pid='plant_id',
+                     function=max_NPV, exclude=None):
+    exclude = exclude if exclude else []
+    with VectorTopo(economic, mode='r') as eco:
+        select_pids = 'SELECT {pid} FROM {tname} GROUP BY {pid};'
+        etab = eco.table
+        exclude.extend((etab.key, eco_pid))
+        ecols = [c for c in etab.columns if c not in exclude]
+        cpids = etab.execute(select_pids.format(pid=eco_pid, tname=etab.name))
+        pids = list(cpids)  # transform the cursor to a list
+        cpids.close()  # close cursor otherwise: dblock error...
+        msgr = get_msgr()
+        with VectorTopo(segment, mode='r') as seg:
+            stab = seg.table
+            scols = set(stab.columns.names())
+            # create the new columns if needed
+            ucols = []
+            msgr.message('Check if column from economic already exists')
+            #import ipdb; ipdb.set_trace()
+            for col in ecols:
+                column = basename + col
+                if column not in scols:
+                    stab.columns.add(column, etab.columns[col])
+                ucols.append(column)
+            for pid, in pids:
+                print('%10s: ' % pid, end='')
+                select_cats = 'SELECT {cat} FROM {tname} WHERE {cpid} LIKE "{pid}"'
+                ecats = list(etab.execute(select_cats.format(cat=etab.key,
+                                                             tname=etab.name,
+                                                             cpid=eco_pid,
+                                                             pid=pid)))
+                print('structures found, ', end='')
+                ec0, ec1 = ecats[0][0], ecats[1][0]
+                l0 = eco.cat(int(ec0), 'lines', layer=1)[0]
+                l1 = eco.cat(int(ec1), 'lines', layer=1)[0]
+                eattrs = function(l0, l1).attrs
+                scats, = list(stab.execute(select_cats.format(cat=stab.key,
+                                                             tname=stab.name,
+                                                             cpid=seg_pid,
+                                                             pid=pid)))
+                if len(scats) != 1:
+                    import ipdb; ipdb.set_trace()
+                print('segment found, ', end='')
+                # TODO: this is not efficient should be done in one step
+                # to avoid to call several time the db update
+                sattr = seg.cat(int(scats[0]), 'lines', layer=1)[0].attrs
+                for ecol, scol in zip(ecols, ucols):
+                    sattr[scol] = str(eattrs[ecol])
+                print('segment updated!')
+            stab.conn.commit()
+            print('Finish')
+
+
+def main(opts, flgs):
+    # check or generate raster map from rules files
+    ecovalues = ['landvalue', 'tributes', 'stumpage', 'rotation', 'age',
+                 'min_exc', 'max_exc']
+    (lan, tri, stu, rot, age,
+     excmin, excmax) = check_raster_or_landuse(opts, ecovalues)
+    upper = opts['upper'] if opts['upper'] else 'tmp_upper'
+    comp = opts['compensation'] if opts['compensation'] else 'tmp_compensation'
+    exc = opts['excavation'] if opts['excavation'] else 'tmp_excavation'
+    vlayer = int(opts['plant_layer'])
+    
+    plant, mset = (opts['plant'].split('@') if '@' in opts['plant'] else (opts['plant'], ''))
+                
+    struct, mset = (opts['struct'].split('@') if '@' in opts['struct'] else (opts['struct'], ''))
+
+    # read common scalar parameters
+    irate = float(opts['interest_rate'])
+    life = float(opts['life'])
+    width = float(opts['width'])
+    depth = float(opts['depth'])
+    slim = float(opts['slope_limit'])
+    overw = overwrite()
+
+    # RASTERS
+    # Start computing the raster map of the costs
+    # Upper value
+    upper_value(upper, stu, lan, rot, age, irate, overw)
+    # Compensation raster costs
+    compensation_cost(comp, lan, tri, upper,
+                      irate, float(opts['gamma_comp']), life, width, overw)
+    # Excavation raster costs
+    excavation_cost(exc, excmin, excmax, opts['slope'],
+                    slim, width, depth, overw)
+    # TODO: extra cost when crossing roads and rivers are missing
+
+    # VECTOR
+    # add columns with costs from rasters
+    # add compensation costs
+    v.rast_stats(map=struct, layer=vlayer, flags='c',
+                 raster=comp, column_prefix='comp_cost', method='sum')
+    # add excavation costs
+    v.rast_stats(map=struct, layer=vlayer, flags='c',
+                 raster=exc, column_prefix='exc_cost', method='sum')
+
+    # add elecro-mechanical costs
+    electromechanical_cost(struct,
+                           power=opts['plant_power_column'],
+                           head=opts['plant_head_column'],
+                           gamma=float(opts['gamma_em']),
+                           alpha=float(opts['alpha_em']),
+                           beta=float(opts['beta_em']),
+                           const=float(opts['const_em']),
+                           vlayer=vlayer, cname='em_cost',
+                           ctype='double precision',
+                           overwrite=overw)
+
+    # add linear cost for pipeline
+    linear_cost(vname=struct, cname='lin_pipe_cost',
+                alpha=float(opts['lc_pipe']), vlayer=vlayer,
+                ctype='double precision',  overwrite=overw)
+
+    # add linear for for electroline
+    get_electro_length(opts)
+    linear_cost(vname=struct, cname='lin_electro_cost',
+                alpha=float(opts['lc_electro']), length='electro_length',
+                vlayer=vlayer, ctype='double precision',  overwrite=overw)
+
+    xcost = "{cname} = {alpha} * {em}"
+    # add power station costs
+    vcolcalc(vname=struct, vlayer=vlayer,
+             ctype='double precision',
+             expr=xcost.format(cname='station_cost', em='em_cost',
+                               alpha=opts['alpha_station']))
+    # add inlet costs
+    vcolcalc(vname=struct, vlayer=vlayer,
+             ctype='double precision', notfinitesubstitute=0.,
+             expr=xcost.format(cname='inlet_cost', em='em_cost',
+                               alpha=opts['alpha_inlet']))
+    # add total inlet costs
+    # TODO: to be check to avoid to count cost more than one time I have moltiplied by 0.5
+    tot = ('tot_cost = (comp_cost_sum + em_cost + '
+           'lin_pipe_cost + lin_electro_cost + '
+           'station_cost + inlet_cost + {grid}*0.5) * '
+           '(1 + {general} + {hindrances})')
+    vcolcalc(vname=struct, vlayer=vlayer,
+             ctype='double precision', notfinitesubstitute=0.,
+             expr=tot.format(grid=opts['grid'], general=opts['general'],
+                             hindrances=opts['hindrances']))
+
+    # TODO: produce a new vector map (output), with the conduct + penstock in
+    # a unique line and sum the costs grouping by intake_id and side
+    # SELECT {key} FROM {tname}
+    #FIXME: intake_id and discharge can have different names
+    group_by(struct, opts['output_struct'],
+             isolate=['intake_id', opts['plant_id_column'],
+                      opts['plant_side_column'], opts['plant_power_column'],
+                      opts['plant_head_column'], 'discharge'],
+             aggregate=['tot_cost', ],
+             function='sum',
+             group_by=['intake_id', opts['plant_side_column']])
+
+    """
+    where these values (3871.2256 and -0.45) are coming from?
+    import numpy as np
+    from scipy import stats
+
+    power= np.array([50., 100., 200., 400., 600., 1000., 5000.])
+    maint = np.array([707., 443., 389., 261., 209., 163., 88.])
+
+    plt.plot(np.log(power), np.log(maint))
+    plt.show()
+    slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(power), np.log(maint))
+    #slope -0.45000431409701719
+    #intercept 8.2613264284076049
+    np.exp(intercept) * power ** slope
+
+    """
+#    maint = "{cname} = {alpha} * {power} ** (1 + {beta}) + {const}"
+    maint = "{cname} = {cost_per_kW} * {power} * {alpha} + {const}"
+    # compute yearly maintenance costs
+    vcolcalc(vname=opts['output_struct'], vlayer=vlayer,
+             ctype='double precision', notfinitesubstitute=0.,
+             expr=maint.format(cname='maintenance',
+                               cost_per_kW=opts['cost_maintenance_per_kw'],
+                               alpha=opts['alpha_maintenance'],
+                               power=opts['plant_power_column'],
+                               beta=opts['beta_maintenance'],
+                               const=opts['const_maintenance']))
+
+    # compute yearly revenues
+    rev = "{cname} = {eta} * {power} * {eprice} * {ophours} * {alpha} + {const}"
+    vcolcalc(vname=opts['output_struct'], vlayer=vlayer,
+             ctype='double precision', notfinitesubstitute=0.,
+             expr=rev.format(cname='revenue',
+                             eta=opts['eta'],
+                             power=opts['plant_power_column'],
+                             eprice=opts['energy_price'],
+                             ophours=opts['operative_hours'],
+                             alpha=opts['alpha_revenue'],
+                             const=opts['const_revenue']))
+
+    # compute the Net Present Value
+    npv = "{cname} = {gamma} * ({revenue} - {maintenance}) - {tot}"
+    gamma_npv = get_gamma_NPV(irate, life)
+    vcolcalc(vname=opts['output_struct'], vlayer=vlayer,
+             ctype='double precision', notfinitesubstitute=0.,
+             expr=npv.format(cname='NPV',
+                             gamma=gamma_npv,
+                             revenue='revenue',
+                             maintenance='maintenance',
+                             tot='tot_cost'))
+
+    economic2segment(economic=opts['output_struct'], segment=plant,
+                         basename=opts['segment_basename'],
+                         eco_layer=1, seg_layer=int(opts['segment_layer']),
+                         eco_pid=opts['plant_id_column'],
+                         seg_pid=opts['segment_column_id'],
+                         function=max_NPV, 
+                         exclude=['intake_id', 'side', 'power', 
+                                  'gross_head', 'discharge'])
+                         
+    vec = VectorTopo(opts['output_struct'])    
+    vec.open('rw')
+    vec.table.columns.add('max_NPV','VARCHAR(3)')
+                                 
+    list_intakeid=list(set(vec.table.execute('SELECT intake_id FROM %s' %vec.table.name).fetchall()))                     
+    
+    for i in range(0,len(list_intakeid)):   
+        vec.rewind()                  
+        list_npv=list(vec.table.execute('SELECT NPV FROM %s WHERE intake_id=%i;' % (vec.table.name, list_intakeid[i][0])).fetchall())
+        npvmax=max(list_npv)[0]
+        for line in vec:
+            if line.attrs['intake_id'] == list_intakeid[i][0]:
+                if line.attrs['NPV'] == npvmax:
+                    line.attrs['max_NPV']='yes'
+                else:
+                    line.attrs['max_NPV']='no'
+
+    vec.table.conn.commit()    
+    vec.close()
+
+
+if __name__ == "__main__":
+    main(*parser())


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r.green.hydro.financial.py
___________________________________________________________________
Added: svn:mime-type
   + text/x-python
Added: svn:eol-style
   + native

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_compmap.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_compmap.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_excmap.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_excmap.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fcompensation.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fcompensation.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fexcavation.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fexcavation.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fgamma.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fgamma.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_flinearcost.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_flinearcost.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fnpv.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fnpv.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fvu.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_fvu.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_input.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_input.png
___________________________________________________________________
Added: svn:mime-type
   + image/png

Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_tableeco.png
===================================================================
(Binary files differ)


Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.financial/r_green_hydro_financial_tableeco.png
___________________________________________________________________
Added: svn:mime-type
   + image/png



More information about the grass-commit mailing list