#!/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@PERMANENT action=prep key_field=cat file_format=unix # v.out.marxan input=test@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())