[GRASS-SVN] r67806 - in grass-addons/grass7/raster/r.futures: . r.futures.potential
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Feb 10 15:19:47 PST 2016
Author: annakrat
Date: 2016-02-10 15:19:47 -0800 (Wed, 10 Feb 2016)
New Revision: 67806
Added:
grass-addons/grass7/raster/r.futures/r.futures.potential/
grass-addons/grass7/raster/r.futures/r.futures.potential/Makefile
grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.html
grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.py
Modified:
grass-addons/grass7/raster/r.futures/Makefile
Log:
r.futures: add Potential submodel
Modified: grass-addons/grass7/raster/r.futures/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/Makefile 2016-02-10 23:11:48 UTC (rev 67805)
+++ grass-addons/grass7/raster/r.futures/Makefile 2016-02-10 23:19:47 UTC (rev 67806)
@@ -6,6 +6,7 @@
r.futures.calib \
r.futures.devpressure \
r.futures.demand \
+ r.futures.potential \
r.futures.pga
include $(MODULE_TOPDIR)/include/Make/Dir.make
Added: grass-addons/grass7/raster/r.futures/r.futures.potential/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.potential/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.potential/Makefile 2016-02-10 23:19:47 UTC (rev 67806)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.futures.potential
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.html (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.html 2016-02-10 23:19:47 UTC (rev 67806)
@@ -0,0 +1,50 @@
+<h2>DESCRIPTION</h2>
+Module <em>r.futures.potential</em> implements POTENTIAL submodel
+as a part of <a href="r.futures.html">FUTURES</a> land change model.
+POTENTIAL is implemented using a set of coefficients that
+relate a selection of site suitability factors to the probability
+of a place becoming developed. This is implemented using
+the parameter table in combination with maps of those site
+suitability factors (mapped predictors).
+
+The coefficients are obtained by conducting multilevel logistic regression in R
+with package <a href="https://cran.r-project.org/web/packages/lme4/index.html">lme4</a>
+where the coefficients may vary by county.
+The best model is selected automatically using <tt>dredge</tt>
+function from package
+<a href="https://cran.r-project.org/web/packages/MuMIn/index.html">MuMIn</a>
+(which has numerous caveats).
+
+<h2>NOTES</h2>
+Note that this module is designed to automate the FUTURES workflow by brute-force
+selection of model, which has numerous caveats.
+
+<h2>EXAMPLES</h2>
+
+<h2>SEE ALSO</h2>
+
+<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.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.sample.category.html">r.sample.category</a></em>
+
+<h2>REFERENCES</h2>
+<ul>
+<li>
+ Meentemeyer, R. K., Tang, W., Dorning, M. A., Vogler, J. B., Cunniffe, N. J., & Shoemaker, D. A. (2013).
+ <a href="http://dx.doi.org/10.1080/00045608.2012.707591">FUTURES: Multilevel Simulations of Emerging
+ Urban-Rural Landscape Structure Using a Stochastic Patch-Growing Algorithm</a>.
+ Annals of the Association of American Geographers, 103(4), 785-807.
+<li>Dorning, M. A., Koch, J., Shoemaker, D. A., & Meentemeyer, R. K. (2015).
+ <a href="http://dx.doi.org/10.1016/j.landurbplan.2014.11.011">Simulating urbanization scenarios reveals
+ tradeoffs between conservation planning strategies</a>.
+ Landscape and Urban Planning, 136, 28-39.</li>
+</ul>
+
+<h2>AUTHOR</h2>
+
+Anna Petrasova, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
+
+<p><i>Last changed: $Date: 2015-02-12 14:32:45 -0500 (Thu, 12 Feb 2015) $</i>
Added: grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.py
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.py (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.py 2016-02-10 23:19:47 UTC (rev 67806)
@@ -0,0 +1,182 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+##############################################################################
+#
+# MODULE: r.futures.potential
+#
+# AUTHOR(S): Anna Petrasova (kratochanna gmail.com)
+#
+# PURPOSE: FUTURES Potential submodel
+#
+# COPYRIGHT: (C) 2016 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: Module for computing development potential as input to r.futures.pga
+#% keyword: raster
+#% keyword: statistics
+#%end
+#%option G_OPT_V_INPUT
+#%end
+#%option G_OPT_F_OUTPUT
+#% description: Output Potential file
+#%end
+#%option G_OPT_DB_COLUMNS
+#% description: Names of attribute columns representing sampled predictors
+#% required: yes
+#%end
+#%option G_OPT_DB_COLUMN
+#% key: developed_column
+#% description: Name of attribute column representing development
+#% required: yes
+#%end
+#%option G_OPT_DB_COLUMN
+#% key: subregions_column
+#% description: Name of attribute column representing subregions
+#% required: yes
+#%end
+#%option
+#% type: integer
+#% key: min_variables
+#% description: Minimum number of predictors considered
+#% required: no
+#% answer: 1
+#% options: 1-20
+#%end
+#%option
+#% type: integer
+#% key: max_variables
+#% description: Maximum number of predictors considered
+#% required: no
+#% options: 1-20
+#%end
+#%flag
+#% key: d
+#% description: Use dredge fuction to find best model
+#%end
+
+import sys
+import atexit
+import subprocess
+import grass.script as gscript
+
+
+rscript = """
+# load required libraries
+for (package in c("MuMIn", "lme4", "optparse")) {
+ if (!(is.element(package, installed.packages()[,1]) ))
+ stop(paste("Package", package, " not found"))
+ }
+suppressPackageStartupMessages(library(MuMIn))
+suppressPackageStartupMessages(library(lme4))
+suppressPackageStartupMessages(library(optparse))
+option_list = list(
+ make_option(c("-i","--input"), action="store", default=NA, type='character', help="input CSV file"),
+ make_option(c("-o","--output"), action="store", default=NA, type='character', help="output CSV file"),
+ make_option(c("-l","--level"), action="store", default=NA, type='character', help="level variable name"),
+ make_option(c("-r","--response"), action="store", default=NA, type='character', help="binary response variable name"),
+ make_option(c("-d","--usedredge"), action="store", default=NA, type='logical', help="use dredge to find best model"),
+ make_option(c("-m","--minimum"), action="store", default=NA, type='integer', help="minimum number of variables for dredge"),
+ make_option(c("-x","--maximum"), action="store", default=NA, type='integer', help="maximum number of variables for dredge")
+)
+
+opt = parse_args(OptionParser(option_list=option_list))
+
+# import data
+input_data = read.csv(opt$input)
+# create global model with all variables
+predictors <- names(input_data)
+predictors <- predictors[predictors != opt$level]
+predictors <- predictors[predictors != opt$response]
+
+interc <- paste("(1|", opt$level, ")")
+fmla <- as.formula(paste(opt$response, " ~ ", paste(c(predictors, interc), collapse= "+")))
+model = glmer(formula=fmla, family = binomial, data=input_data, na.action = "na.fail")
+
+if(opt$usedredge) {
+ #create all possible models, always include county as the level
+ select.model <- dredge(model, evaluate=TRUE, rank="AIC", fixed=~(1|opt$level), m.lim=c(opt$minimum, opt$maximum), trace=FALSE)
+
+ # save the best model
+ model.best <- get.models(select.model, 1)
+ model = glmer(formula(model.best[[1]]), family = binomial, data=input_data, na.action = "na.fail")
+}
+print(summary(model))
+coefs <- as.data.frame(coef(model)[[1]])
+write.table(cbind(rownames(coefs), coefs), opt$output, row.names=FALSE, sep="\t")
+"""
+
+TMP_CSV = None
+TMP_POT = None
+TMP_RSCRIPT = None
+
+
+def cleanup():
+ gscript.try_remove(TMP_CSV)
+ gscript.try_remove(TMP_POT)
+ gscript.try_remove(TMP_RSCRIPT)
+
+
+def main():
+ vinput = options['input']
+ columns = options['columns'].split(',')
+ binary = options['developed_column']
+ level = options['subregions_column']
+ minim = int(options['min_variables'])
+ dredge = flags['d']
+ if options['max_variables']:
+ maxv = (options['max_variables'])
+ else:
+ maxv = len(columns) - 2
+ if dredge and minim > maxv:
+ gscript.fatal(_("Minimum number of predictor variables is larger than maximum number"))
+
+ global TMP_CSV, TMP_RSCRIPT, TMP_POT
+ TMP_CSV = gscript.tempfile(create=False) + '.csv'
+ TMP_RSCRIPT = gscript.tempfile()
+ with open(TMP_RSCRIPT, 'w') as f:
+ f.write(rscript)
+ TMP_POT = gscript.tempfile(create=False) + '_potential.csv'
+
+ columns += [binary, level]
+ gscript.run_command('v.db.select', map=vinput, columns=columns, separator='comma', file=TMP_CSV)
+
+ if dredge:
+ gscript.info(_("Running automatic model selection ..."))
+ else:
+ gscript.info(_("Computing model..."))
+ p = subprocess.Popen(['Rscript', TMP_RSCRIPT, '-i', TMP_CSV, '-l', level, '-r', binary, '-m', str(minim), '-x', str(maxv), '-o', TMP_POT, '-d', 'TRUE' if dredge else 'FALSE'],
+ stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+ stdout, stderr = p.communicate()
+ print stderr
+ if p.returncode != 0:
+ print stderr
+ gscript.fatal(_("Running R script failed, check messages above"))
+
+ gscript.info(_("Best model summary:"))
+ gscript.info("-------------------------")
+ print stdout
+
+ with open(TMP_POT, 'r') as fin, open(options['output'], 'w') as fout:
+ i = 0
+ for line in fin.readlines():
+ row = line.strip().split('\t')
+ row = [each.strip('"') for each in row]
+ if i == 0:
+ row[0] = "ID"
+ row[1] = "Intercept"
+ fout.write(' '.join(row))
+ fout.write('\n')
+ i += 1
+
+
+if __name__ == "__main__":
+ options, flags = gscript.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
More information about the grass-commit
mailing list