[GRASS-SVN] r44441 - in grass-addons/vector: . v.out.marxan

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Nov 26 17:21:19 EST 2010


Author: tsw
Date: 2010-11-26 14:21:18 -0800 (Fri, 26 Nov 2010)
New Revision: 44441

Added:
   grass-addons/vector/v.out.marxan/
   grass-addons/vector/v.out.marxan/Makefile
   grass-addons/vector/v.out.marxan/description.html
   grass-addons/vector/v.out.marxan/v.out.marxan
Log:
New script for creating Marxan input files

Added: grass-addons/vector/v.out.marxan/Makefile
===================================================================
--- grass-addons/vector/v.out.marxan/Makefile	                        (rev 0)
+++ grass-addons/vector/v.out.marxan/Makefile	2010-11-26 22:21:18 UTC (rev 44441)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.out.marxan
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

Added: grass-addons/vector/v.out.marxan/description.html
===================================================================
--- grass-addons/vector/v.out.marxan/description.html	                        (rev 0)
+++ grass-addons/vector/v.out.marxan/description.html	2010-11-26 22:21:18 UTC (rev 44441)
@@ -0,0 +1,47 @@
+<H2>DESCRIPTION</H2>
+
+<EM>v.out.marxan</EM> is a python script prepare a vecotr planning grid file and create marxan input files. 
+Marxan output can be read back into GRASS for display with <EM>v.in.marxan</EM>. 
+
+<H2>NOTES</H2>
+
+Use of the DOS output format will create a second file with '.dos' appended to the output file name.  
+Consistent specification of the Planning Unit Key field is necessary in all Marxan scripts to ensure a working set of ouptut files. 
+All GRASS Marxan modules require Python, PostgreSQL 8.x or higher and are intended for use with Marxan version 1.8.10. 
+
+One of the options for output is the scaling of boundary and species variables to ranges between 0 and 1. Since the 
+objective function of Marxan is sum of the planning unit cost plus the sum of the boundary cost * the boundary length 
+modifier plus the sum of species penalty factors, it is desirable that all of these be kept within the same order 
+of magnitude. The scaling function for data outputs is designed to address this.
+
+The second potentially useful feature is in the calculation of the boundary length. A Calculated option is provided
+that multiplies the actual length in map units by the specified attribute of the adjacent planning unit. In this case
+the user needs to specify whether the greater or lesser value will be applied in all cases where they differ. This calculated
+value is then scaled to range from 0 to 1. This allows for using attributes attached to each planning unit to 
+define the relative merit of joining them together. This might be particularly useful in cases where certain habitats
+have synergetic effects and others do not work well together. 
+
+These brief notes do not in any constitute a outline of Marxan's capabilities or provide guidance to appropriate use. 
+For this please familiarize yourself with the user manual and good practices handbook found on the Marxan home page.
+
+<h2>SEE ALSO</h2>
+
+<em><a HREF="v.out.marxan.html">v.out.marxan</a>,
+<a HREF="v.mkhexgrid">v.mkhexgrid</a>,
+</em>
+<p><a HREF="http://www.uq.edu.au/marxan/">Marxan Home Page</a></p>
+
+<H2>AUTHOR</H2>
+Trevor Wiens
+
+<H2>REFERENCES</H2>
+
+<p><a href="http://www.uq.edu.au/marxan/docs/marxan_manual_1_8_2.pdf" target="_blank">Ball, I.R., and H.P. Possingham, 2000. MARXAN (V1.8.2): Marine Reserve Design Using Spatially Explicit Annealing, a Manual.</a></p>
+
+<p><a href="http://www.uq.edu.au/marxan/docs/Marxan_User_Manual_2008.pdf" target="_blank">Game, E.T. and H.S. Grantham, 2008. Marxan User Manual: For Marxan version 1.8.10. University of Queensland, St. Lucia, Queensland, Australia, and Pacific Marine Analysis and Research Association, Vancouver, British Columbia, Canada.</a></p>
+
+<p><a href="http://www.uq.edu.au/marxan/docs/Possingham_etal_2000_Mathematical.pdf" target="_blank">Possingham, H. P., I. R. Ball and S. Andelman (2000) Mathematical methods for identifying representative reserve networks. In: S. Ferson and M. Burgman (eds) Quantitative methods for conservation biology. Springer-Verlag, New York, pp. 291-305.</a></p>
+
+<p><a href="http://www.uq.edu.au/marxan/docs/Marxan%20Good%20Practices%20Handbook%20v2%202010.pdf" target="_blank">Ardron, J. H.P. Possingham and C.J. Klein (Eds.),Version 2, 2010. Marxan good practices handbook. University of Queensland, St. Lucia, Queensland, Australia, and Pacific Marine Analysis and Research Association, Vancouver, British Columbia, Canada.</a> </p>
+
+<p><i>Last changed: $Date: 2010-11-24 $</i>

Added: grass-addons/vector/v.out.marxan/v.out.marxan
===================================================================
--- grass-addons/vector/v.out.marxan/v.out.marxan	                        (rev 0)
+++ grass-addons/vector/v.out.marxan/v.out.marxan	2010-11-26 22:21:18 UTC (rev 44441)
@@ -0,0 +1,737 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE:	v.out.marxan for GRASS 6.4 (2010-07-14)
+#               modified TSW (2010-11-22)
+#		
+# AUTHOR(S):	Trevor Wiens 
+#
+# PURPOSE: Creates the following tabs do the following
+#       required tab - choose what to do or ouptut to generate
+#       boundary tab - bound.dat
+#       features tab - spec.dat and puvspr2.dat
+#       planning units tab  - pu.dat
+#
+# REQUIREMENTS: PostgreSQL 8.x or above, Marxan 1.8.10
+#
+# COPYRIGHT:	(C) Trevor Wiens, 2010 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.
+#
+# EXAMPLES:
+#       v.out.marxan input=test at PERMANENT action=prep key_field=cat file_format=unix
+#       v.out.marxan input=test at PERMANENT action=boundary key_field=cat file_format=unix bound_cfield=cat
+#############################################################################
+
+#%Module
+#%  description: Prepares vector file and generates output for use in Marxan.
+#%  keywords: vector
+#%  keywords: marxan
+#%End
+#%option
+#%  key: input
+#%  type: string
+#%  gisprompt: old,vector,vector
+#%  label: Planning Unit Vector
+#%  description: Name of planning unit vector
+#%  required: yes
+#%end
+#%option
+#%  key: action
+#%  type: string
+#%  options: prep,boundary,features,planning,all
+#%  label: Action
+#%  description: Prepare planning unit vector, export boundary file, export feature files, export planning file, or all. 
+#%  multiple: no
+#%  answer: prep
+#%  required: yes
+#%end
+#%option
+#%  key: key_field
+#%  type: string
+#%  label: Planning Unit Key or ID Field
+#%  description: Name of key field for planning unit
+#%  answer: cat
+#%  required: yes
+#%end
+#%option
+#%  key:file_format
+#%  type: string
+#%  options: unix,dos
+#%  label: Output File Format
+#%  description: Format of output file.
+#%  multiple: no
+#%  answer: unix
+#%  required: yes
+#%  guisection: Main
+#%end
+
+#%option
+#%  key: bound_out
+#%  type: string
+#%  gisprompt: new_file,file,output
+#%  answer: bound.dat
+#%  label: Boundary File Name
+#%  description: Name of boundary output file. The Marxan default is bound.dat.
+#%  required : no
+#%  guisection: Boundary
+#%end
+#%option
+#%  key:bound_length
+#%  type: string
+#%  options: fixed,m,km,scaled,field,calculated
+#%  label: Boundary Length Options
+#%  description: Fixed length, metres, kilometres, scaled from 0 to 1, use value of adjacent PUs, or scaled value of length * value of adjacent PUs.
+#%  multiple: no
+#%  answer: scaled
+#%  required: no
+#%  guisection: Boundary
+#%end
+#%option
+#%  key: bound_cfield
+#%  type: string
+#%  answer: 
+#%  label: PU Field Name
+#%  description: Field value of adjacent PU's used to adjust boundary length
+#%  required : no
+#%  guisection: Boundary
+#%end
+#%option
+#%  key: bound_caction
+#%  type: string
+#%  options: greatest,least
+#%  answer: greatest
+#%  label: Adjacent PU Value Action
+#%  description: Select which field value for adjacent planning units to use to modify the boundary lengths.
+#%  required : no
+#%  guisection: Boundary
+#%end
+
+#%option
+#%  key: feature_spec_out
+#%  type: string
+#%  gisprompt: new_file,file,output
+#%  answer: spec.dat
+#%  label: Conservation Feature File Name
+#%  description: Name of conservation feature or species output file. The Marxan default is spec.dat.
+#%  required : no
+#%  guisection: Features
+#%end
+#%option
+#%  key: feature_vs_out
+#%  type: string
+#%  gisprompt: new_file,file,output
+#%  answer: puvspr2.dat
+#%  label: PU vs Feature File Name
+#%  description: Name of planning unit vs conservation feature (species) output file name. The Marxan default is puvspr2.dat.
+#%  required : no
+#%  guisection: Features
+#%end
+#%option
+#%  key: feature_out_format
+#%  type: string
+#%  options: raw,scaled
+#%  label: Conservation Factor Format
+#%  description: Format of conservation factor as either raw values or scaled from 0 to 1.
+#%  multiple: no
+#%  answer: scaled
+#%  required: no
+#%  guisection: Features
+#%end
+#%option
+#%  key: feature_fields
+#%  type: string
+#%  answer: 
+#%  multiple: yes
+#%  label: Conservation Factor Field Names
+#%  description: List conservation factor field names separated by commas.
+#%  required: no
+#%  guisection: Features
+#%end
+#%option
+#%  key: feature_names
+#%  type: string
+#%  answer: 
+#%  multiple: yes
+#%  label: Conservation Factor Names 
+#%  description: List of descriptive names for the conservation factors in the same order as above separated by commas. 
+#%  guisection: Features
+#%end
+#%option
+#%  key: feature_penalties
+#%  type: string
+#%  answer: 
+#%  multiple: yes
+#%  label: Conservation Factor Penalty Factors
+#%  description: List the penalty factors for each conservation feature in the same order above separated by commas.
+#%  required: no
+#%  guisection: Features
+#%end
+#%option
+#%  key: feature_targets
+#%  type: string
+#%  answer: 
+#%  multiple: yes
+#%  label: Conservation Factor Targets
+#%  description: List conservatoion targets for each feature in the same order as above separated by commas. 
+#%  required: no
+#%  guisection: Features
+#%end
+
+#%option
+#%  key: plan_out
+#%  type: string
+#%  gisprompt: new_file,file,output
+#%  answer: pu.dat
+#%  label: Planning Unit File Name
+#%  description: Name of planning unit output file. The Marxan default is pu.dat.
+#%  required : no
+#%  guisection: Planning Units
+#%end
+#%option
+#%  key: plan_cost
+#%  type: string
+#%  options: none,fixed,field,field_scaled
+#%  label: Planning Unit Cost
+#%  description: Cost of planning unit as none, a fixed value, raw field value, or field scaled from 0 to 1.
+#%  multiple: no
+#%  answer: field_scaled
+#%  required: no
+#%  guisection: Planning Units
+#%end
+#%option
+#%  key: plan_cost_field
+#%  type: string
+#%  label: Planning Unit Cost Field Name
+#%  description: Name of field used to calculate planning unit costs.
+#%  required: no
+#%  guisection: Planning Units
+#%end
+#%option
+#%  key: plan_status
+#%  type: string
+#%  options: no_assumptions,initial_inclusion,field
+#%  label: Planning Unit Status
+#%  description: Initial conditions for planning units set as no assumptions, initial inclusion, or set by field content.
+#%  multiple: no
+#%  answer: no_assumptions
+#%  required: no
+#%  guisection: Planning Units
+#%end
+#%option
+#%  key: plan_status_field
+#%  type: string
+#%  label: Planning Unit Status Field Name
+#%  description: Name of field used to set intial status of planning unit fields.
+#%  required: no
+#%  guisection: Planning Units
+#%end
+
+import sys
+import os
+import subprocess
+import time
+from grass.script import core as grass
+from grass.script import vector as gvect
+from grass.script import db as gdb
+
+#
+# Common Functions
+# 
+
+def Unix2DOS(unix_text):
+    try:
+        test = unix_text.index('\r')
+        dos_text = unix_text
+    except:
+        dos_text = unix_text.replace('\n','\r\n')
+    return(dos_text)
+    
+def writeDOSfile(sourcefile):
+    dosfname = sourcefile +'.dos'
+    tmpfs = file(sourcefile, 'r')
+    stext = tmpfs.read()
+    tmpfs.close()
+    otext = Unix2DOS(stext)
+    tmpfd = file(dosfname,'w')
+    tmpfd.write(otext)
+    tmpfd.close()
+
+def clean_up(tempname):
+    grass.run_command('g.remove', vect=tempname,quiet=True)
+
+#
+# Vector File Preparation
+#
+
+#
+# prepVector - checks to see if pu_status is present, if so assumes file is properly prepared.
+#              if pu_status is not present, adds this along with xloc, yloc, mxn_best, and mxn_freq fields
+#              these are all set to zero and then xloc and yloc are updated with pu centroid locations
+#
+# modified TSW - 22 November 2010: Replaced use of numeric for double and direct queries with v.db.addcol to 
+#                                  facilitate backend database flexibility.
+#       
+
+def prepVector():
+    table_name=gvect.vector_db(options['input'])[1]['table']
+    # use pipes to capture errors
+    errpipe = subprocess.PIPE
+    outpipe = subprocess.PIPE
+    # test if preparation already executed
+    querytext = 'select pu_status from %s' % table_name
+    #grass.message(querytext)
+    try:
+        r = gdb.db_select(table=table_name, sql=querytext, flags='-q', stdout = outpipe, stderr = errpipe)
+    except:
+        # a problem means we need to prepare the vector
+        pass
+    else:
+        # no problem means abort
+        grass.message('Vector %s already prepared' % table_name)
+        return(0)
+    # create columns
+    new_cols = 'pu_status int, xloc double precision, yloc double precision, mxn_best int, mxn_freq int' 
+    querytext = 'update %s set pu_status = 0, xloc = 0.0, yloc = 0.0, mxn_best = 0, mxn_freq = 0' % table_name
+    try:
+        grass.run_command('v.db.addcol', map='%s' % options['input'], columns=new_cols)
+        r = grass.write_command('db.execute', stdin = querytext, stdout = outpipe, stderr = errpipe)
+        if r <> 0:
+            raise
+    except:
+        msgtext="An error occured when altering the input table. \n" + \
+            "Please review query and table structure for possible conflicts \n"  + querytext
+        grass.message(msgtext, flag='e')
+        return(-1)
+    # notify success
+    grass.message('Columns created.')
+    #cmd = "v.to.db map=%s type=centroid layer=1 qlayer=1 option=coor col=xloc,yloc" % options['input']
+    try:
+        r = grass.run_command('v.to.db', map='%s' % options['input'], type='centroid', \
+            layer=1, qlayer=1, option='coor', col='xloc,yloc', quiet=True)
+    except:
+        grass.message("An error occurred updating field values.", flag='e')
+        return(-1)
+    # notify success
+    if r == 0:
+        grass.message("Vector File Preparation Complete")
+    return(0)
+     
+
+#
+# Boundary File Export
+#
+
+#
+# exportBoundary - Checks to make sure all the needed fields are completed before passing on work to 
+#                  calculate_adjancency_length function.
+#
+#
+
+def exportBoundary():
+    if (options['bound_length'] == 'calculated' or options['bound_length'] == 'field') and options['bound_cfield'] == '':
+        grass.message('Field required if calculated or field values are selected')
+        return(-1)
+    core_table_name = gvect.vector_db(options['input'])[1]['table']
+    # grass.message("Core Table: %s" % core_table_name)
+    tempname= 'tempvect%d' % os.getpid()
+    testval = calculate_adjacency_length(tempname, core_table_name)
+    if testval == -1:
+        return(-1)
+    bound_table_name = gvect.vector_db(tempname)[2]['table']
+    # grass.message("Boundary Table: %s" % bound_table_name)    
+    grass.message("Writing boundary file")
+    stat = write_boundary_file(core_table_name,bound_table_name)
+    clean_up(tempname)
+    if stat == -1:
+        return(-1)
+    else:
+        return(0)
+
+#
+# calculate_adjacency_length - Uses simple series of grass commands to calculate adjacency and length 
+#                              of all planning units.
+#
+# parameters:
+#               tempname - name of temporary file
+#               corename - name of main vector attribute table
+#
+
+def calculate_adjacency_length(tempname, corename):
+    grass.message("Creating temporary file")
+    #
+    # create temporary boundary layer file
+    #
+    grass.run_command('v.category', input=options['input'], output=tempname, layer='2', \
+        type='boundary', option='add', quiet=True)
+    # determine database schema
+    info = grass.read_command('v.db.connect', flags = 'g', map=tempname, key='cat', layer=1, fs='|')
+    outline = info.split('|')[1]
+    fullname = outline.split('.')
+    if len(fullname) > 1:
+        schemaname = outline.split('.')[0]
+    else:
+        # if no schema assume public (postgresql default)
+        schemaname = 'public'
+    tablename = schemaname + '.'  + tempname + '_2'
+    #
+    # add attribute table to hold boundary adjacency and length information
+    #
+    # create table
+    new_cols = 'left_cat int, right_cat int, link_cat int, s_length double precision, left_cval double precision, '
+    new_cols += 'right_cval double precision, cval double precision'
+    grass.run_command('v.db.addtable', map=tempname, table=tablename , layer='2', columns=new_cols, quiet=True)
+    # update the information
+    grass.message("Updating information")
+    grass.run_command('v.to.db', map=tempname, layer='2', option='sides', columns='left_cat,right_cat', quiet=True)
+    if options['bound_length'] == 'km':
+        grass.run_command('v.to.db', map=tempname, layer='2', option='length', units='k', columns='s_length', quiet=True)
+    else:
+        grass.run_command('v.to.db', map=tempname, layer='2', option='length', units='me', columns='s_length', quiet=True)
+    # NOTE - the following SQL code is written for PostgreSQL and will probably fail with other database backends
+    if options['bound_length'] == 'calculated' or options['bound_length'] == 'field':
+        # update the left value
+        sqltext1 = """
+            update %s
+            set left_cval = b.%s
+            from %s b
+            where %s.left_cat = b.cat;""" % (tablename, options['bound_cfield'], corename, tablename)
+        # update the right value
+        sqltext2 = """
+            update %s
+            set right_cval = b.%s
+            from %s b
+            where %s.right_cat = b.cat;""" % (tablename, options['bound_cfield'], corename, tablename)
+        # set cval to zero
+        sqltext3 = "update %s set cval = 0;" % tablename
+        # conditionally set cval based on greatest or least option
+        sqltext4 = """
+            update %s
+            set cval = %s(left_cval, right_cval)
+            where left_cat <> -1 and right_cat <> -1;""" % (tablename, options['bound_caction'])
+        try:
+            errpipe = subprocess.PIPE
+            outpipe = subprocess.PIPE
+            errnum = 1
+            grass.write_command('db.execute', stdin=sqltext1, stdout = outpipe, stderr = errpipe)
+            errnum = 2
+            grass.write_command('db.execute', stdin=sqltext2, stdout = outpipe, stderr = errpipe)
+            errnum = 3
+            grass.write_command('db.execute', stdin=sqltext3, stdout = outpipe, stderr = errpipe)
+            errnum = 4
+            grass.write_command('db.execute', stdin=sqltext4, stdout = outpipe, stderr = errpipe)
+        except:
+            grass.message("Errnum: %s Field (%s) is not valid. Please ensure that field exists in source table" % (errnum, options['bound_cfield']))
+            grass.message("sql: %s " % sqltext1)
+            return(-1)
+
+#
+# write_boundary_file - this finishes the calculations and writes the results to a text file
+#                      
+# parameters:
+#               core - name of main vector attribute table
+#               bound - name of boundary attribute table
+# 
+# modified TSW - 22 November 2010: Adjusted to ensure that border units have their boundary lengths included to 
+#                                  prevent articially low values and thus selection of borders over other PUs
+#
+
+def write_boundary_file(core, bound):
+    # select the correct query text
+    if options['bound_length'] == 'm' or options['bound_length'] == 'km':
+	sqltext = """select a.left_cat as id1,
+		(case when a.right_cat = -1 then a.left_cat else a.right_cat end) as id2,
+		a.s_length as boundary
+		from (
+		select distinct left_cat, right_cat, sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat
+		order by left_cat, right_cat) a
+	""" %  (bound)
+    elif options['bound_length'] == 'scaled':
+	sqltext = """select a.left_cat as id1,
+		(case when a.right_cat = -1 then a.left_cat else a.right_cat end) as id2,
+		a.s_length / b.maxlen as boundary
+		from 
+		(select distinct left_cat, right_cat, sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat
+		order by left_cat, right_cat) a,
+                (select max(q.s_length) as maxlen
+                from (select distinct left_cat, right_cat, sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat
+		order by left_cat, right_cat) q ) b
+	""" %  (bound, bound)
+    elif options['bound_length'] == 'field':
+	sqltext = """select a.left_cat as id1,
+		(case when a.right_cat = -1 then a.left_cat else a.right_cat end) as id2,
+		a.cval as boundary
+		from (
+		select distinct left_cat, right_cat, cval sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat, cval
+		order by left_cat, right_cat, cval) a
+	""" %  (bound)
+    elif options['bound_length'] == 'calculated':
+	sqltext = """select a.left_cat as id1,
+		(case when a.right_cat = -1 then a.left_cat else a.right_cat end) as id2,
+		(a.cval * a.s_length) / b.maxlen as boundary
+		from (
+		select distinct left_cat, right_cat, cval sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat, cval
+		order by left_cat, right_cat, cval) a,
+                (select max(z.boundary) as maxlen 
+                from
+                (select q.left_cat as id1,
+		(case when q.right_cat = -1 then q.left_cat else q.right_cat end) as id2,
+		q.cval * q.s_length as boundary
+		from (
+		select distinct left_cat, right_cat, cval sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat, cval
+		order by left_cat, right_cat, cval) q ) z ) b
+	""" %  (bound,bound)
+    else:
+        # fixed value
+	sqltext = """select a.left_cat as id1,
+		(case when a.right_cat = -1 then a.left_cat else a.right_cat end) as id2,
+		1 as boundary
+		from (
+		select distinct left_cat, right_cat, sum(s_length) as s_length
+		from %s
+		where not(left_cat = -1 and right_cat = -1)
+		group by left_cat, right_cat
+		order by left_cat, right_cat) a
+	""" %  (bound)
+    # test query and write file
+    try:
+        grass.run_command('db.select', flags='t', sql=sqltext, quiet=True)
+        tmpf = file(options['bound_out'], 'w')
+        grass.run_command('db.select', sql=sqltext, fs='\t',stdout = tmpf)
+        tmpf.close()
+    except:
+        grass.message("Boundary export query failed not execute. \nPlease check the field names are correct. ", flag='e')
+        return(-1)
+    # convert to dos if needed
+    if options['file_format'] == 'dos':
+        writeDOSfile(options['bound_out'])
+    return(0)
+    
+#
+# Feature Files Export
+#
+
+#
+# exportFeatures - a wrapper function that exports both puvspr2.dat and spec.dat
+#
+#
+
+def exportFeatures():
+    table_name=gvect.vector_db(options['input'])[1]['table']
+    field_list = options['feature_fields'].split(',')
+    stat = write_vs_file(table_name, field_list)
+    if stat == -1:
+        return(stat)
+    write_species_file(table_name, field_list)
+    return(0)
+    
+#
+# write_vs_file - writes puvspr2.dat
+#
+# parameters:
+#               table_name - name of vector attribute table
+#               field_list - name of vector fields to be exported as planning unit attributes
+#
+# NOTE: This function depends on PostgreSQL because of use of grass interface can not over-ride formatting
+#       of large values into scientific notation which fails in Marxan.
+#    
+    
+def write_vs_file(table_name, field_list):
+    # create query for each field and 
+    # union them together
+    utext = ''
+    for i in range(0,len(field_list)):
+        if options['feature_out_format'] == 'raw':
+            qtext = """
+                select %d as species, a.%s as pu, to_char(a.%s::numeric, 
+                'FM999999999999999999999999999999.9999999999999999999') as amount
+                from %s a
+            """ % (i, options['key_field'],field_list[i], table_name)
+        elif options['feature_out_format'] == 'scaled':
+            qtext = """
+                select %d as species, a.%s as pu, to_char(a.%s::numeric / b.spec_max::numeric,
+                'FM999999999999999999999999999999.9999999999999999999') as amount
+                from %s a,
+                (select max(%s) as spec_max
+                from %s) b
+            """ % (i, options['key_field'],field_list[i], table_name, field_list[i], table_name)
+        if utext == '':
+            utext = qtext
+        else:
+            utext = utext + ' union ' + qtext
+    # embed union text into larger query
+    sqltext = """select m.species, m.pu, m.amount
+        from (%s) m
+        where m.amount::numeric > 0
+        order by m.species, m.pu""" % utext
+    #print sqltext
+    # test query and create output file
+    #grass.message(sqltext)
+    try:
+        grass.run_command('db.select', flags='t', sql=sqltext)
+        tmpf = file(options['feature_vs_out'], 'w')
+        grass.run_command('db.select', sql=sqltext, fs='\t', stdout=tmpf)
+        tmpf.close()
+    except:
+        grass.message("Feature export query failed not execute. \nPlease check the field names are correct.", flag='e')
+        return(-1)
+    # create dos version if necessary
+    if options['file_format'] == 'dos':
+        writeDOSfile(options['feature_vs_out'])
+
+#
+# write_species_file - writes spec.dat aka feat.dat
+#
+# parameters:
+#               table_name - vector attribute table name
+#               field_list - list of attribute fields to be exported as conservation features or species
+#
+
+def write_species_file(table_name, field_list):
+    name_list = options['feature_names'].split(',')
+    spf_list = options['feature_penalties'].split(',')
+    target_list = options['feature_targets'].split(',')
+    tmpf = file(options['feature_spec_out'], 'w')
+    tmpf.write("id\ttarget\tspf\tname\n")
+    for i in range(0,len(field_list)):
+        tmpf.write("%d\t%s\t%s\t%s\n" % (i,target_list[i],spf_list[i],name_list[i]))
+    tmpf.close()
+    # create dos version if necessary
+    if options['file_format'] == 'dos':
+        writeDOSfile(options['feature_spec_out'])
+
+#
+# Planning Unit Export
+#
+
+#
+# exportPlanningUnits - exports planning units file (pu.dat)
+#
+# NOTE: This function depends on PostgreSQL casting functions
+#
+
+def exportPlanningUnits():
+    # get table name containing vector attributes
+    table_name=gvect.vector_db(options['input'])[1]['table']
+    if options['plan_cost'] == 'none':
+        # no cost layer
+        if options['plan_status'] == 'no_assumptions':
+            sqltext = "select %s as id, 0 as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], table_name, options['key_field'])
+        elif options['plan_status'] == 'initial_inclusion':
+            sqltext = "select %s as id, 1 as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], table_name, options['key_field'])
+        else:
+            sqltext = "select %s as id, %s as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], options['plan_status_field'], table_name, options['key_field'])
+    elif options['plan_cost'] == 'fixed':
+        # all units the same value
+        if options['plan_status'] == 'no_assumptions':
+            sqltext = "select %s as id, 1 as cost, 0 as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], table_name, options['key_field'])
+        elif options['plan_status'] == 'initial_inclusion':
+            sqltext ="select %s as id, 1 as cost, 1 as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], table_name, options['key_field'])
+        else:
+            sqltext = "select %s as id, 1 as cost, %s as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], options['plan_status_field'], table_name, options['key_field'])
+    elif options['plan_cost'] == 'field':
+        # planning unit costs defined by named field
+        if options['plan_status'] == 'no_assumptions':
+            sqltext = "select %s as id, %s as cost, 0 as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], options['plan_cost_field'], table_name, options['key_field'])
+        elif options['plan_status'] == 'initial_inclusion':
+            sqltext = "select %s as id, %s as cost, 1 as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], options['plan_cost_field'], table_name, options['key_field'])
+        else:
+            sqltext = "select cat as id, %s as cost, %s as status, xloc, yloc from %s order by %s" % \
+            (options['key_field'], options['plan_cost_field'], options['plan_status_field'], table_name, options['key_field'])
+    elif options['plan_cost'] == 'field_scaled':
+        # planning unit costs defined by named field and scaled from 0 to 1
+        if options['plan_status'] == 'no_assumptions':
+            sqltext = """select a.%s as id, cast(a.%s as numeric) / cast(b.scale as numeric) as cost, 
+                0 as status, xloc, yloc from %s a, 
+                (select max(%s) as scale from %s) b
+                order by %s
+                """ % (options['key_field'], options['plan_cost_field'], table_name, \
+                options['plan_cost_field'], table_name, options['key_field'])
+        elif options['plan_status'] == 'initial_inclusion':
+            sqltext = """select a.%s as id, cast(a.%s as numeric) / cast(b.scale as numeric) as cost, 
+                1 as status, xloc, yloc from %s a, 
+                (select max(%s) as scale from %s) b
+                order by %s
+                """ % (options['key_field'], options['plan_cost_field'], table_name, \
+                options['plan_cost_field'], table_name, options['key_field'])
+        else:
+            sqltext = """select a.%s as id, cast(a.%s as numeric) / cast(b.scale as numeric) as cost, 
+                %s as status, xloc, yloc from %s a, 
+                (select max(%s) as scale from %s) b
+                order by %s
+                """ % (options['key_field'], options['plan_cost_field'], options['plan_status_field'], table_name, \
+                options['plan_cost_field'], table_name, options['key_field'])
+    #grass.message(sqltext)
+    try:
+        grass.run_command('db.select', flags='t', sql=sqltext)
+        tmpf = file(options['plan_out'], 'w')
+        grass.run_command('db.select', sql=sqltext, fs='\t', stdout=tmpf)
+        tmpf.close()
+    except:
+        grass.message("Planning export query failed not execute. \nPlease check the field names are correct.", flag='e')
+        return(-1)
+    
+    if options['file_format'] == 'dos':
+        writeDOSfile(options['output'])
+    return(0)
+
+#
+# Main function
+#
+    
+def main():
+    starttime = time.localtime()    
+    if options['action'] == 'prep' or options['action'] == 'all':
+        status = prepVector()
+        if status == -1:
+            return(-1)
+    if options['action'] == 'boundary' or options['action'] == 'all':
+        status = exportBoundary()
+        if status == -1:
+            return(-1)
+    if options['action'] == 'features' or options['action'] == 'all':
+        status = exportFeatures()
+        if status == -1:
+            return(-1)
+    if options['action'] == 'planning' or  options['action'] == 'all':
+        status = exportPlanningUnits()
+        if status == -1:
+            return(-1)
+    endtime = time.localtime()
+    msgtext ='Done!\nStarted: '+ time.strftime("%Y.%m.%d %H:%M:%S",starttime) + \
+    '\nFinished:'+ time.strftime("%Y.%m.%d %H:%M:%S",endtime)
+    grass.message(msgtext)
+    return(0)
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    sys.exit(main())


Property changes on: grass-addons/vector/v.out.marxan/v.out.marxan
___________________________________________________________________
Added: svn:executable
   + *



More information about the grass-commit mailing list