[GRASS-SVN] r60804 - grass-addons/grass6/raster/r.floorsim
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 12 11:25:02 PDT 2014
Author: isaacullah
Date: 2014-06-12 11:25:02 -0700 (Thu, 12 Jun 2014)
New Revision: 60804
Added:
grass-addons/grass6/raster/r.floorsim/r.floorsim
Log:
A script to simulate post-depositional "site formation processes" in archaeological house-floor contexts.
Added: grass-addons/grass6/raster/r.floorsim/r.floorsim
===================================================================
--- grass-addons/grass6/raster/r.floorsim/r.floorsim (rev 0)
+++ grass-addons/grass6/raster/r.floorsim/r.floorsim 2014-06-12 18:25:02 UTC (rev 60804)
@@ -0,0 +1,211 @@
+#!/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()
+
Property changes on: grass-addons/grass6/raster/r.floorsim/r.floorsim
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list