[GRASS-SVN] r71316 - grass-addons/grass7/raster/r.futures/r.futures.potential
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jul 26 20:15:52 PDT 2017
Author: annakrat
Date: 2017-07-26 20:15:52 -0700 (Wed, 26 Jul 2017)
New Revision: 71316
Modified:
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
Log:
r.futures.potential: allow using only one subregion
Modified: 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 2017-07-25 07:36:30 UTC (rev 71315)
+++ grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.html 2017-07-27 03:15:52 UTC (rev 71316)
@@ -26,7 +26,8 @@
The header contains the names of the predictor maps and the first column
contains the identifiers of the subregions. The order of columns is important,
the second column represents intercept, the third development pressure and
-then the predictors.
+then the predictors. Therefore the development pressure column must be specified
+as the first column in option <b>columns</b>.
<pre>
ID Intercept devpressure_0_5 slope road_dens_perc forest_smooth_perc ...
@@ -41,6 +42,8 @@
<h2>NOTES</h2>
Note that this module is designed to automate the FUTURES workflow by brute-force
selection of model, which has numerous caveats.
+<p>
+In case there is only one subregion, R function glm is used instead of glmer.
<h2>EXAMPLES</h2>
Modified: 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 2017-07-25 07:36:30 UTC (rev 71315)
+++ grass-addons/grass7/raster/r.futures/r.futures.potential/r.futures.potential.py 2017-07-27 03:15:52 UTC (rev 71316)
@@ -92,23 +92,41 @@
input_data = read.csv(opt$input)
# create global model with all variables
predictors <- names(input_data)
-predictors <- predictors[predictors != opt$level]
+if (!is.na(opt$level)) {
+ 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 (is.na(opt$level)) {
+ fmla <- as.formula(paste(opt$response, " ~ ", paste(c(predictors), collapse= "+")))
+ model = glm(formula=fmla, family = binomial, data=input_data, na.action = "na.fail")
+} else {
+ 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)
-
+ if (is.na(opt$level)) {
+ select.model <- dredge(model, evaluate=TRUE, rank="AIC", m.lim=c(opt$minimum, opt$maximum), trace=FALSE)
+ } else {
+ 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")
+ if (is.na(opt$level)) {
+ model = glm(formula(model.best[[1]]), family = binomial, data=input_data, na.action = "na.fail")
+ } else {
+ 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]])
+if (is.na(opt$level)) {
+ coefs <- t(as.data.frame(coef(model)))
+} else {
+ coefs <- as.data.frame(coef(model)[[1]])
+}
write.table(cbind(rownames(coefs), coefs), opt$output, row.names=FALSE, sep="\t")
"""
@@ -133,18 +151,25 @@
if options['max_variables']:
maxv = (options['max_variables'])
else:
- maxv = len(columns) - 2
+ maxv = len(columns)
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()
+ include_level = True
+ distinct = gscript.read_command('v.db.select', flags='c', map=vinput,
+ columns="distinct {l}".format(l=level)).strip()
+ if len(distinct.splitlines()) <= 1:
+ include_level = False
+ single_level = distinct.splitlines()[0]
with open(TMP_RSCRIPT, 'w') as f:
f.write(rscript)
TMP_POT = gscript.tempfile(create=False) + '_potential.csv'
-
- columns += [binary, level]
+ columns += [binary]
+ if include_level:
+ columns += [level]
where = "{c} IS NOT NULL".format(c=columns[0])
for c in columns[1:]:
where += " AND {c} IS NOT NULL".format(c=c)
@@ -154,8 +179,12 @@
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)
+
+ cmd = ['Rscript', TMP_RSCRIPT, '-i', TMP_CSV, '-r', binary,
+ '-m', str(minim), '-x', str(maxv), '-o', TMP_POT, '-d', 'TRUE' if dredge else 'FALSE']
+ if include_level:
+ cmd += [ '-l', level]
+ p = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
print stderr
if p.returncode != 0:
@@ -174,6 +203,8 @@
if i == 0:
row[0] = "ID"
row[1] = "Intercept"
+ if i == 1 and not include_level:
+ row[0] = single_level
fout.write('\t'.join(row))
fout.write('\n')
i += 1
@@ -183,3 +214,4 @@
options, flags = gscript.parser()
atexit.register(cleanup)
sys.exit(main())
+
More information about the grass-commit
mailing list