#!/usr/bin/python ############################################################################ # # MODULE: r.floorsim.py # AUTHOR(S): iullah # COPYRIGHT: (C) 2014, Isaac Ullah # # description: Simulation of artifacts deposition on a housefloor, with site formation disturbance, sampling, and re-randomization. # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # #############################################################################/ #%Module #% description: Simulation of artifacts deposition on a housefloor, with site formation disturbance, sampling, and re-randomization. #%End #%option #% key: prefix #% type: string #% description: Prefix for all output maps and files. #% answer: floorsim_ #% required : yes #%END #%option #% key: areas #% type: string #% gisprompt: old,cell,raster #% description: Map of activity areas in which to deposit artifacts #% required : yes #%END #%option #% key: density #% type: string #% description: Resolution of "artifact" deposition (in meters, density will be 1 artifact per cell at this cell resolution) #% answer: 0.01 #% required : yes #%END #%option #% key: disturbance #% type: string #% description: Maximum amount of movement (meters) of artifacts due to site formation processes (actual distance will vary along a normal distribution around the origin) #% answer: 0.5 #% required : yes #%END #%option #% key: s_points #% type: string #% gisprompt: old,vector,vector #% description: Vector points map of sampling locations (table will be updated, must be in current mapset) #% required : yes #%END #%option #% key: s_size #% type: string #% description: Size of the samplng squares (diameter of samples) centered on the sampling points (meters) #% answer: 0.1 #% required : yes #%END #%option #% key: runs #% type: string #% description: Number of re-randomized runs (montecarlo) to do #% answer: 100 #% required : yes #%END #%Flag #% key: c #% description: Only output stats (do not save any maps) #%end import os import sys import random import grass.script as grass def main(): prefix = options['prefix'] areas = options['areas'] density = options['density'] disturbance = options['disturbance'] s_points = options['s_points'] s_size = options['s_size'] runs= int(options['runs']) flag_c = flags['c'] env = grass.gisenv() #set the region ot match the arctivity areas map grass.run_command('g.region', flags = 'a', rast = areas) #test to make sure sampling points map is in current mapset, otherwise return a fatal exit. if s_points.split('@')[0] not in grass.list_grouped('vect')[env['MAPSET']]: grass.fatal('Run Terminated!\n\n' + s_points + ' MUST be in current mapset.\n\nPlease copy ' + s_points + ' into current mapset and try again!') #open a text file to write stats to, and write headers stats_name = prefix + 'stats_file.csv' statsfile = file(stats_name, 'w') statsfile.write('Run_number,Covariance,R2,1st_quartile_pos_deviation,Median_pos_deviation,3rd_quartile_pos_deviation,1st_quartile_neg_deviation,Median_neg_deviation,3rd_quartile_neg_deviation\n') #Start a loop to do the iterative randomized runs for x in range(runs): number = str(x+1) #set the tracking number for this run grass.message('Run # %s:\n-----------------' % number) seednum = random.randrange(1000) #generate a random seed number for this run #set file names for this run init_pts = prefix + 'initial_points_' + number dist_pts = prefix + 'disturbed_points_' + number dens_map = prefix + 'real_density_map_' + number smo_map = prefix + 'smoothed_real_density_map' + number col ='run_' + number #Test to see if colum name exists. If it does, return error message with fatal exit for item in grass.db_describe(s_points.split('@')[0])['cols']: if col in item: grass.fatal('Sorry, table ' + s_points.split('@')[0] + ' already contains a column named ' + col+ '\nPlease delete this column, or create a new copy of the sampling points map that does not have this column, and try again.' ) interp_map = prefix + 'interpolated_density_map_' + number dif_map = prefix + 'difference_map_' + number pos_dev = prefix + 'positive_deviations_map_' + number neg_dev = prefix + 'negative_deviations_map_' + number #start the process for this run grass.message("initializing run...") #set the resolution to the density for initial depostion of "artifacts" grass.run_command('g.region', flags = 'a', res = density) grass.message('depositing artifacts...') #deposit "artifacts" evenly (one/cell) at given density across the activity areas grass.run_command('r.random', quiet = True, input = areas, n = '100%', vector_output = init_pts) grass.message('site formation disturbance...') # disturb the locations of the "artifacts" to simulate site formation proccesses (this is the core engine for differences between runs) grass.run_command('v.perturb', quiet = True, input = init_pts, output = dist_pts, distribution = 'normal', parameters = '0,' + disturbance, seed = seednum) grass.message('creating map of actual artifact density...') # set the region to the sampling diameter, so that we collect density info at a comparable resolution. grass.run_command('g.region', flags = 'a', res = s_size) #create a map of the real density at resolution od sampling units grass.run_command('v.neighbors', quiet = True, input = dist_pts, output = dens_map, method = 'count', size = s_size) #change null values to 0 so that the stats are right grass.run_command('r.null', quiet = True, map = dens_map, null = '0') #Since the interpolated maps are by their nature smoothed, we need to smooth the raw real density maps to be able to compare them with the interpolated maps. grass.run_command('r.neighbors', quiet = True, input = dens_map, output = smo_map, size = 5, method = 'average') grass.message('computing artfiact density at sampling points...') #create a new column in the samplingpoints file to update with the density values from this run colname = 'run_' + number + ' integer' grass.run_command('v.db.addcolumn', quiet = True, map = s_points, columns = colname) #write density values at sampling points to the table grass.run_command('v.what.rast', quiet =True, vector = s_points, raster = dens_map, column =col) grass.message('creating interpolated map from sampled data...') #make interpolated density map from sampled data from this run grass.run_command('v.surf.rst', quiet = True, input = s_points, zcolumn = col, elev = interp_map, tension = '60', smooth = '2', segmax = '600') grass.message('calculating statistics for this run...') #grab the coviance between real and terpolated density maps cv = grass.pipe_command('r.covar', flags = 'r', map = smo_map + ',' + interp_map) cv_dict = {} for line in cv.stdout: [key, val] = line.strip().split() cv_dict[key] = val cv.wait() #figure otu the r-squared value for linear regression between real and interpolated density maps reg = grass.pipe_command('r.regression.line', flags = 'gs', map1 = smo_map, map2 = interp_map) reg_dict = {} for line in reg.stdout: [key, val] = line.strip().split('=') reg_dict[key] = val reg.wait() r2 = str(pow(float(reg_dict['R']), 2)) #create map of differences between real density map and interpolated density map grass.mapcalc("${out} = ${rast1} - ${rast2}", out = dif_map, rast1 = smo_map, rast2 = interp_map) #make a map of all the positive values (for stats) grass.mapcalc("${out} = if(${rast} > 0, ${rast}, null())", out = pos_dev, rast = dif_map) #make a map of all the negative values (for stats) grass.mapcalc("${out} = if(${rast} < 0, ${rast}, null())", out = neg_dev, rast = dif_map) #grab the stats for the positive deviations (areas that were overpredicted by the interpolation) pd = grass.pipe_command('r.univar', flags = 'ge', map = pos_dev) pd_stats = {} for line in pd.stdout: [key, val] = line.strip().split('=') pd_stats[key] = val pd.wait() #grab the stats for the negative deviations (areas that were underpredicted by the interpolation) nd = grass.pipe_command('r.univar', flags = 'ge', map = neg_dev) nd_stats = {} for line in nd.stdout: [key, val] = line.strip().split('=') nd_stats[key] = val nd.wait() grass.message('writing stats to text file...') #write this run's stats to the stats file statsfile.write(number + ',' + cv_dict['1.000000'] + ',' + r2 + ',' + pd_stats['first_quartile'] + ',' + pd_stats['median'] + ',' + pd_stats['third_quartile'] + ',' + nd_stats['first_quartile'] + ',' + nd_stats['median'] + ',' + nd_stats['third_quartile'] + '\n') #if we keep maps, we color them, otherwise, we delete them! if flag_c: grass.run_command('g.mremove', quiet = True, flags = 'f', rast = prefix + '*' + number, vect = prefix + '*' + number) else: grass.run_command('r.colors', quiet = True, map = interp_map, raster = dens_map) grass.run_command('r.colors', quiet = True, map = dif_map, color = 'differences') grass.run_command('r.colors', quiet = True, map = pos_dev, color = 'bgyr') grass.run_command('r.colors', quiet = True, flags = 'n', map = neg_dev, color = 'bgyr') grass.message('-----------------') #The loop is done, now close the stats file statsfile.close() grass.message('-------------------\nALL RUNS COMPLETE\n-------------------') return 0 if __name__ == "__main__": options, flags = grass.parser() main()