/
r.floorsim
executable file
·211 lines (200 loc) · 10.3 KB
/
r.floorsim
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
#!/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()