[GRASS-SVN] r66582 - in grass-addons/grass7/raster/r.futures: . r.futures.calib r.futures.demand r.futures.devpressure r.futures.pga
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Oct 22 18:50:11 PDT 2015
Author: annakrat
Date: 2015-10-22 18:50:11 -0700 (Thu, 22 Oct 2015)
New Revision: 66582
Added:
grass-addons/grass7/raster/r.futures/r.futures.demand/
grass-addons/grass7/raster/r.futures/r.futures.demand/Makefile
grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.html
grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.py
grass-addons/grass7/raster/r.futures/r.futures.demand/r_futures_demand_plot.png
Modified:
grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html
grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html
grass-addons/grass7/raster/r.futures/r.futures.html
grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html
Log:
r.futures: add demand submodule
Modified: grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html 2015-10-22 22:21:46 UTC (rev 66581)
+++ grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html 2015-10-23 01:50:11 UTC (rev 66582)
@@ -92,6 +92,7 @@
<a href="r.futures.html">FUTURES</a>,
<em><a href="r.futures.pga.html">r.futures.pga</a></em>,
<em><a href="r.futures.devpressure.html">r.futures.devpressure</a></em>,
+<em><a href="r.futures.demand.html">r.futures.demand</a></em>,
<em><a href="r.sample.category.html">r.sample.category</a></em>
<h2>REFERENCES</h2>
Added: grass-addons/grass7/raster/r.futures/r.futures.demand/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.demand/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.demand/Makefile 2015-10-23 01:50:11 UTC (rev 66582)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.futures.demand
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.futures/r.futures.demand/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.html (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.html 2015-10-23 01:50:11 UTC (rev 66582)
@@ -0,0 +1,102 @@
+<h2>DESCRIPTION</h2>
+The <em>r.futures.demand</em> module of FUTURES
+determines the quantity of expected land changed.
+It creates a demand table as the number of cells to be converted
+at each time step for each subregion based on the relation
+between the population and developed land in the past years.
+
+<p>
+The input accepts multiple (at least 2)
+rasters of developed (category 1) and undeveloped areas (category 0)
+from different years, ordered by time.
+For these years, user has to provide the population numbers for each subregion
+in parameter <b>observed_population</b> as a CSV file.
+The format is as follows. First column is time
+(matching the time of rasters used in parameter <b>development</b>) and
+first row is the category of the subregion.
+<pre>
+year 1 2 ...
+1985 19860 10980 ...
+1995 20760 12660 ...
+2005 21070 13090 ...
+2015 22000 13940 ...
+</pre>
+<p>
+The same table is needed for projected population
+(parameter <b>projected_population</b>).
+The categories of the input raster <b>subregions</b> must
+match the identifiers of subregions in files given in <b>observed_population</b>
+and <b>projected_population</b>.
+Parameter <b>simulation_times</b> is a comma separated list
+of times for which the demand will be computed. The first time should
+be the time of the developed/undeveloped raster used
+in <em><a href="r.futures.pga">r.futures.pga</a></em>
+as a starting point for simulation. There is an easy way to create
+such list using Python:
+<pre>
+','.join([str(i) for i in range(2015, 2031)])
+</pre>
+or Bash:
+<pre>
+seq -s, 2015 2030
+</pre>
+
+<p>
+The <b>method</b> parameter allows to choose the type of relation
+between population and developed area. The available methods
+include linear, logarithmic and exponential relation.
+If more than one method is checked, the best relation is selected
+based on RMSE. Recommended methods are logarithmic or linear.
+
+<p>
+The format of the output <b>demand</b> table is:
+<pre>
+Years_to_simulate: 24
+year 1 2 ...
+2016 2900 243 ...
+2017 2790 240 ...
+2018 2660 232 ...
+2019 2345 221 ...
+2020 2154 215 ...
+...
+</pre>
+where each value represents the number of new developed cells in each step.
+In case the demand values would be negative (in case of population decrease
+or if the relation is inversely proportional) the values are turned into zeros,
+since FUTURES does not simulate change from developed to undeveloped sites.
+By changing the "Years_to_simulate" number, user can control the simulation.
+
+<p>
+An optional output <b>plot</b> is a plot of the relations for each subregion.
+It allows to more effectively assess the relation suitable for each subregion.
+
+<center>
+<img src="r_futures_demand_plot.png">
+<p>
+Figure: Example of logarithmic relation between population and developed area
+(generated with option <b>plot</b>).
+</center>
+
+<h2>NOTES</h2>
+<em>r.futures.demand</em> computes the relation
+between population and developed area using simple
+regression. It therefore serves as starting point
+for demand computation which can be further improved.
+
+<h2>EXAMPLES</h2>
+<div class="code"><pre>
+r.futures.demand devel=ash_lc_1976,ash_lc_1985,ash_lc_1995,ash_lc_2006 \
+sub=index2_in observ_pop=ashvile_observed_population.csv \
+proj_pop=ashvile_projected_population.csv simul=`seq -s, 2007 2030` sep=comma \
+demand=demand_sub1.txt plot=regression_lin.pdf method=linear
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.futures.pga.html">r.futures.pga</a></em>
+
+<h2>AUTHOR</h2>
+
+Anna Petrasova, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.py
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.py (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.py 2015-10-23 01:50:11 UTC (rev 66582)
@@ -0,0 +1,240 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+##############################################################################
+#
+# MODULE: r.futures.demand
+#
+# AUTHOR(S): Anna Petrasova (kratochanna gmail.com)
+#
+# PURPOSE: create demand table for FUTURES
+#
+# COPYRIGHT: (C) 2015 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (version 2). Read the file COPYING that comes with GRASS
+# for details.
+#
+##############################################################################
+
+#%module
+#% description: Script for creating demand table which determines the quantity of land change expected.
+#% keyword: raster
+#% keyword: demand
+#%end
+#%option G_OPT_R_INPUTS
+#% key: development
+#% description: Names of input binary raster maps representing development
+#% guisection: Input maps
+#%end
+#%option G_OPT_R_INPUT
+#% key: subregions
+#% description: Raster map of subregions
+#% guisection: Input maps
+#%end
+#%option G_OPT_F_INPUT
+#% key: observed_population
+#% description: CSV file with observed population in subregions at certain times
+#% guisection: Input population
+#%end
+#%option G_OPT_F_INPUT
+#% key: projected_population
+#% description: CSV file with projected population in subregions at certain times
+#% guisection: Input population
+#%end
+#%option
+#% type: integer
+#% key: simulation_times
+#% multiple: yes
+#% required: yes
+#% description: For which times demand is projected
+#% guisection: Output
+#%end
+#%option
+#% type: string
+#% key: method
+#% multiple: yes
+#% required: yes
+#% description: Relationship between developed cells (dependent) and population (explanatory)
+#% options: linear, logarithmic, exponential
+#% descriptions:linear;y = A + Bx;logarithmic;y = A + Bln(x);exponential;y = Ae^(BX)
+#% answer: linear,logarithmic
+#% guisection: Optional
+#%end
+#%option G_OPT_F_OUTPUT
+#% key: plot
+#% required: no
+#% description: Save plotted relationship between developed cells and population into a file
+#% guisection: Output
+#%end
+#%option G_OPT_F_OUTPUT
+#% key: demand
+#% description: Output CSV file with demand (times as rows, regions as columns)
+#% guisection: Output
+#%end
+#%option G_OPT_F_SEP
+#% label: Separator used in input CSV files
+#% guisection: Input population
+#%end
+
+
+import sys
+import numpy as np
+
+import grass.script.core as gcore
+import grass.script.utils as gutils
+
+
+def main():
+ developments = options['development'].split(',')
+ observed_popul_file = options['observed_population']
+ projected_popul_file = options['projected_population']
+ sep = gutils.separator(options['separator'])
+ subregions = options['subregions']
+ methods = options['method'].split(',')
+ plot = options['plot']
+ simulation_times = [float(each) for each in options['simulation_times'].split(',')]
+
+ observed_popul = np.genfromtxt(observed_popul_file, dtype=float, delimiter=sep, names=True)
+ projected_popul = np.genfromtxt(projected_popul_file, dtype=float, delimiter=sep, names=True)
+ year_col = observed_popul.dtype.names[0]
+ observed_times = observed_popul[year_col]
+ year_col = projected_popul.dtype.names[0]
+ projected_times = projected_popul[year_col]
+
+ if len(developments) != len(observed_times):
+ gcore.fatal(_("Number of development raster maps doesn't not correspond to the number of observed times"))
+
+ # gather developed cells in subregions
+ gcore.info(_("Computing number of developed cells..."))
+ table_developed = {}
+ subregionIds = set()
+ for i in range(len(observed_times)):
+ gcore.percent(i, len(observed_times), 1)
+ data = gcore.read_command('r.univar', flags='gt', zones=subregions, map=developments[i])
+ for line in data.splitlines():
+ stats = line.split('|')
+ if stats[0] == 'zone':
+ continue
+ subregionId, developed_cells = int(stats[0]), int(stats[12])
+ subregionIds.add(subregionId)
+ if i == 0:
+ table_developed[subregionId] = []
+ table_developed[subregionId].append(developed_cells)
+ gcore.percent(1, 1, 1)
+ subregionIds = sorted(list(subregionIds))
+ # linear interpolation between population points
+ population_for_simulated_times = {}
+ for subregionId in table_developed.keys():
+ population_for_simulated_times[subregionId] = np.interp(x=simulation_times,
+ xp=np.append(observed_times, projected_times),
+ fp=np.append(observed_popul[str(subregionId)],
+ projected_popul[str(subregionId)]))
+ # regression
+ demand = {}
+ i = 0
+ if plot:
+ import matplotlib.pyplot as plt
+ n_plots = np.ceil(np.sqrt(len(subregionIds)))
+ fig = plt.figure(figsize=(5 * n_plots, 5 * n_plots))
+
+ for subregionId in subregionIds:
+ i += 1
+ rmse = dict()
+ predicted = dict()
+ simulated = dict()
+ coeff = dict()
+ for method in methods:
+ # observed population points for subregion
+ reg_pop = observed_popul[str(subregionId)]
+ if method == 'logarithmic':
+ reg_pop = np.log(reg_pop)
+ if method == 'exponential':
+ y = np.log(table_developed[subregionId])
+ else:
+ y = table_developed[subregionId]
+ A = np.vstack((reg_pop, np.ones(len(reg_pop)))).T
+ m, c = np.linalg.lstsq(A, y)[0] # y = mx + c
+ coeff[method] = m, c
+ simulated[method] = np.array(population_for_simulated_times[subregionId])
+ if method == 'logarithmic':
+ predicted[method] = np.log(simulated[method]) * m + c
+ r = (reg_pop * m + c) - table_developed[subregionId]
+ elif method == 'exponential':
+ predicted[method] = np.exp(m * simulated[method] + c)
+ r = np.exp(m * reg_pop + c) - table_developed[subregionId]
+ else: # linear
+ predicted[method] = simulated[method] * m + c
+ r = (reg_pop * m + c) - table_developed[subregionId]
+ # RMSE
+ rmse[method] = np.sqrt((np.sum(r * r) / (len(reg_pop) - 2)))
+
+ method = min(rmse, key=rmse.get)
+ # write demand
+ demand[subregionId] = predicted[method]
+ demand[subregionId] = np.diff(demand[subregionId])
+ if np.any(demand[subregionId] < 0):
+ gcore.warning(_("Subregion {sub} has negative numbers"
+ " of newly developed cells, changing to zero".format(sub=subregionId)))
+ demand[subregionId][demand[subregionId] < 0] = 0
+ if m < 0:
+ demand[subregionId].fill(0)
+ gcore.warning(_("For subregion {sub} population and development are inversely proportional,"
+ "will result in zero demand".format(sub=subregionId)))
+
+ # draw
+ if plot:
+ ax = fig.add_subplot(n_plots, n_plots, i)
+ ax.set_title(str(subregionId) + ', RMSE: ' + str(rmse[method]))
+ ax.set_xlabel('population')
+ ax.set_ylabel('developed cells')
+ # plot known points
+ x = np.array(observed_popul[str(subregionId)])
+ y = np.array(table_developed[subregionId])
+ ax.plot(x, y, marker='o', linestyle='', markersize=8)
+ # plot predicted curve
+ x_pred = np.linspace(np.min(x),
+ np.max(np.array(population_for_simulated_times[subregionId])), 10)
+ m, c = coeff[method]
+ if method == 'linear':
+ line = x_pred * m + c
+ label = "$y = {c:.3f} + {m:.3f} x$".format(m=m, c=c)
+ elif method == 'logarithmic':
+ line = np.log(x_pred) * m + c
+ label = "$y = {c:.3f} + {m:.3f} \ln(x)$".format(m=m, c=c)
+ else:
+ line = np.exp(x_pred * m + c)
+ label = "$y = {c:.3f} e^{{{m:.3f}x}}$".format(m=m, c=np.exp(c))
+ ax.plot(x_pred, line, label=label)
+ ax.plot(simulated[method], predicted[method], linestyle='', marker='o', markerfacecolor='None')
+ plt.legend(loc=0)
+ labels = ax.get_xticklabels()
+ plt.setp(labels, rotation=30)
+ if plot:
+ plt.tight_layout()
+ fig.savefig(plot)
+
+ # write demand
+ with open(options['demand'], 'w') as f:
+ f.write('Years_to_simulate: {sim}\n'.format(sim=len(simulation_times)))
+ f.write('\t'.join(observed_popul.dtype.names))
+ f.write('\n')
+ i = 0
+ for time in simulation_times[1:]:
+ f.write(str(int(time)))
+ f.write('\t')
+ # put 0 where there are more counties but are not in region
+ for sub in subregionIds:
+ if sub not in subregionIds:
+ f.write('0')
+ else:
+ f.write(str(int(demand[sub][i])))
+ if sub != subregionIds[-1]:
+ f.write('\t')
+ f.write('\n')
+ i += 1
+
+
+if __name__ == "__main__":
+ options, flags = gcore.parser()
+ sys.exit(main())
Property changes on: grass-addons/grass7/raster/r.futures/r.futures.demand/r.futures.demand.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.futures/r.futures.demand/r_futures_demand_plot.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.futures/r.futures.demand/r_futures_demand_plot.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Modified: grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html 2015-10-22 22:21:46 UTC (rev 66581)
+++ grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html 2015-10-23 01:50:11 UTC (rev 66582)
@@ -59,6 +59,7 @@
<a href="r.futures.html">FUTURES</a>,
<em><a href="r.futures.pga.html">r.futures.pga</a></em>,
<em><a href="r.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.futures.demand.html">r.futures.demand</a></em>,
<em><a href="r.sample.category.html">r.sample.category</a></em>
<h2>REFERENCES</h2>
Modified: grass-addons/grass7/raster/r.futures/r.futures.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.html 2015-10-22 22:21:46 UTC (rev 66581)
+++ grass-addons/grass7/raster/r.futures/r.futures.html 2015-10-23 01:50:11 UTC (rev 66582)
@@ -20,8 +20,8 @@
on user's preferences and data availability.
Land area conversion over time can be derived for the USA, e.g.
from National Land Cover Dataset.
- A possible implementation of the DEMAND submodel will be available as module
- <em>r.futures.demand</em>.</dd>
+ A possible implementation of the DEMAND submodel is available as module
+ <em><a href="r.futures.demand.html">r.futures.demand</em>.</dd>
<dt><em>POTENTIAL</em></dt>
<dd>The POTENTIAL submodel uses site suitability modeling approaches
@@ -127,6 +127,7 @@
<em><a href="r.futures.pga.html">r.futures.pga</a></em>,
<em><a href="r.futures.devpressure.html">r.futures.devpressure</a></em>,
<em><a href="r.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.futures.demand.html">r.futures.demand</a></em>,
<em><a href="r.sample.category.html">r.sample.category</a></em>
Modified: grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html 2015-10-22 22:21:46 UTC (rev 66581)
+++ grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html 2015-10-23 01:50:11 UTC (rev 66582)
@@ -138,6 +138,7 @@
<a href="r.futures.html">FUTURES</a>,
<em><a href="r.futures.devpressure.html">r.futures.devpressure</a></em>,
<em><a href="r.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.futures.demand.html">r.futures.demand</a></em>,
<em><a href="r.sample.category.html">r.sample.category</a></em>
More information about the grass-commit
mailing list