· 8 years ago · Dec 25, 2017, 02:14 AM
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
61 61
62 62
63 63
64 64
65 65
66 66
67 67
68 68
69 69
70 70
71 71
72 72
73 73
74 74
75 75
76 76
77 77
78 78
79 79
80 80
81 81
82 82
83 83
84 84
85 85
86 86
87 87
88 88
89 89
90 90
91 91
92 92
93 93
94 94
95 95
96 96
97 97
98 98
99 99
100100
101101
102102
103103
104104
105105
106106
107107
108108
109109
110110
111111
112112
113113
114114
115115
116116
117117
118118
119119
120120
121121
122122
123123
124124
125125
126126
127127
128128
129129
130130
131131
132132
133133
134134
135135
136136
137137
138138
139139
140140
141141
142142
143143
144144
145145
146146
147147
148148
149149
150150
151151
152152
153153
154
155
156
157#!/usr/bin/env python
158#cvaSamplesGen
159#Create random points, calculate viewsheds per random sample and patch the results into a file.
160#Viewshed sizes in (n) number of points and (n) number of samples Generator ver. 1.3
161#Useful for statistical analysis
162#The raster file used to generate random points could be different from the one used by r.cva
163#Israel Hinojosa Balino (UCL)
164#This script use r.cva by Mark Lake (UCL)
165#This script needs GRASS 6.4 and Python to run
166#Run this script from command line in the GRASS Environment
167#Reference: Conolly J. and Lake, M. (2006) Geographical Information Systems in Archaeology.
168#Cambridge University Press: New York
169
170import os
171import sys
172grass_install_tree=os.getenv('GISBASE')
173sys.path.append(grass_install_tree+os.sep+'etc'+os.sep+'python')
174import grass.script as grass
175
176class bcolors:
177 HEADER = '\033[7m'
178 OKBLUE = '\033[94m'
179 OKGREEN = '\033[92m'
180 WARNING = '\033[2m'
181 CONFIRM = '\033[1m'
182 FAIL = '\033[94m'
183 ENDC = '\033[0m'
184
185 def disable(self):
186 self.HEADER = ''
187 self.OKBLUE = ''
188 self.OKGREEN = ''
189 self.WARNING = ''
190 self.FAIL = ''
191 self.ENDC = ''
192
193print ("\n" + bcolors.HEADER + bcolors.FAIL + "cvaSamplesGen 1.2" + bcolors.ENDC)
194print (bcolors.WARNING + "Viewshed sizes in (n) number of points and (n) number of samples generator" + bcolors.ENDC)
195print (bcolors.WARNING + "Create viewsheds at random points in a definite number of samples" + "\n"
196 "and patch the results into a file" + bcolors.ENDC + "\n" + "Please write the information required" + "\n")
197print (bcolors.OKGREEN + "Ideally, the extent of the raster used to produce the random points should consider" + "\n"
198 "the radius of the maximum viewing distance to avoid the edge effect" + bcolors.ENDC + "\n")
199
200raster=raw_input("Name of initial raster to calculate random points: ")
201
202#check twice if name of DEM or DTM (raster file) exists, otherwise return error
203filex=grass.find_file(raster, element='cell')
204if not filex['fullname'] != '':
205 raster=raw_input("Raster file with the name"+" "+"<"+ raster +">"+" "+"does not exist. Another try? : ")
206filex=grass.find_file(raster, element='cell')
207if not filex['fullname'] != '':
208 raster=raw_input("Seriously"+" "+"<"+ raster +">"+" "+"does not exist. Please choose another name: ")
209filex=grass.find_file(raster, element='cell')
210if not filex['fullname'] != '':
211 grass.fatal(_("You have to check your files before using this script. Bye!!!"))
212
213print (bcolors.OKBLUE + "\n" + "\"When the distance between a given viewpoint and the edge of the map region" + "\n"
214 "is less than that radius it follows that the viewshed may be artificially" + "\n"
215 "truncated, thus invalidating comparison with the viewsheds of other viewpoints" + "\n"
216 "that were further from the edge of the map\" (Conolly & Lake, 229)" + bcolors.ENDC)
217
218rasterCVA=raw_input("\n" + "Name of raster to calculate the viewshed sizes points with r.cva.: ")
219print (bcolors.OKGREEN + "The extent of this raster should consider a buffer zone" + "\n"
220 "of the same width as the maximum viewing distance" + bcolors.ENDC)
221
222#check twice if name of DEM or DTM (raster file) exists, otherwise return error
223filecva=grass.find_file(rasterCVA, element='cell')
224if not filecva['fullname'] != '':
225 rasterCVA=raw_input("Raster file with the name"+" "+"<"+ rasterCVA +">"+" "+"does not exist. Another try? : ")
226filecva=grass.find_file(rasterCVA, element='cell')
227if not filecva['fullname'] != '':
228 rasterCVA=raw_input("Seriously"+" "+"<"+ rasterCVA +">"+" "+"does not exist. Please choose another name: ")
229filecva=grass.find_file(rasterCVA, element='cell')
230if not filecva['fullname'] != '':
231 grass.fatal(_("You have to check your files before using this script. Bye!!!"))
232
233file=raw_input("\n" + "Name of the output vector file containing samples: ")
234
235#check twice if name of previous output file exists, otherwise return error
236filey=grass.find_file(file, element='vector')
237if filey['fullname'] != '':
238 file=raw_input("A file named"+" "+"<"+ file +">"+" "+"exists. Please choose another name: ")
239filey=grass.find_file(file, element='vector')
240if filey['fullname'] != '':
241 file=raw_input("Seriously"+" "+"<"+ file +">"+" "+"exist. Please choose another name: ")
242filey=grass.find_file(file, element='vector')
243if filey['fullname'] != '':
244 grass.fatal(_("You have to choose another name. Bye!"))
245
246#customization
247max=int(raw_input("\n" + "Number of samples: ") )
248maxR=int(raw_input("Number of random points: ") )
249maxDist=int(raw_input("Maximum distance (in metres) from viewing point (refer to r.cva): ") )
250
251print "\n" + bcolors.HEADER + "Final Check" + bcolors.ENDC
252print "Input raster to generate random points: " + bcolors.CONFIRM + raster + bcolors.ENDC
253print "Input raster to generate viewsheds: " + bcolors.CONFIRM + rasterCVA + bcolors.ENDC
254print "Output file name (vector): " + bcolors.CONFIRM + file + bcolors.ENDC
255print "Number of samples: " + bcolors.CONFIRM, max, bcolors.ENDC
256print "Number of random points per sample " + bcolors.CONFIRM, maxR, bcolors.ENDC
257print "Maximum viewing distance " + bcolors.CONFIRM, maxDist, bcolors.ENDC
258
259yes = raw_input("Continue (yes/no): ")
260if yes == "yes":
261 print ("Perfect!")
262else:
263 grass.fatal(_("Sorry. Bye!"))
264
265#Create empty file, attach a table and add columns with the proper names and types
266grass.run_command('v.edit', map=file, type='line', tool='create')
267grass.run_command('v.db.addtable', map=file)
268grass.run_command('v.db.addcol', map=file, columns="value DOUBLE PRECISION,vsize INT,sample INT")
269
270#generate random points and calculate viewsheds
271i=1
272while i <= max:
273 print i
274 namerandom="randomV" +str(i)
275 nameRandomCVA= "randomCVA"+str(i)
276 sampleNo=i
277 print 'Generating'+' '+'random'+' '+'numbers'+' '+'for'+' '+ namerandom
278 i=i+1
279 grass.run_command('r.random', overwrite='True', input=raster, vector_output=namerandom, n=maxR)
280 print 'Calculating'+' '+'viewsheds'+' '+'for'+' '+ nameRandomCVA
281 grass.run_command('r.cva', 'o', input=rasterCVA, output=nameRandomCVA, sites=namerandom,
282 obs_elev='1.7', target_elev='0.0', max_dist=maxDist, seed='1',
283 sample='10.0',type='sites',curvc='0.0' )
284
285#Add column to attach Viewshed sizes and perform a spatial query to extract the values from raster to vector
286 print 'Spatial'+' '+'query'+' '+'for'+' '+ nameRandomCVA
287 grass.run_command('v.db.addcol', map=namerandom, layer='1', columns='vsize INT')
288 grass.run_command('v.what.rast', vector=namerandom, raster=nameRandomCVA, layer='1', column='vsize')
289
290#Add another column to attach the sample number
291 print 'Adding'+' '+'column'+' '+'for'+' '+'sample'+' '+'numbers'+' '+'in'+' '+ namerandom
292 grass.run_command('v.db.addcol', map=namerandom, layer='1', columns='sample INT')
293 grass.run_command('v.db.update', map=namerandom, layer='1', column='sample', value=sampleNo)
294
295#Append sample to pre-existing samples
296 print 'Patching'+' '+'original'+' '+'vector'+' '+'with'+' '+ namerandom
297 grass.run_command('v.patch', 'ae',overwrite='True', input=namerandom, output=file)
298
299#Clean dataset from temporal files
300 print 'Deleting'+' '+'supporting'+' '+'files'+' '+':'+' '+ namerandom +' '+'and'+' '+ nameRandomCVA
301 grass.run_command('g.remove', vect=namerandom, rast=nameRandomCVA)
302 print 'DONE!'
303
304#export to ascii with coordinates
305 print 'Generating ASCII file with headers and coordinates'
306 grass.run_command('v.out.ascii', input=file, output='samples.txt', dp=10, columns='vsize,sample')
307 print 'DONE!'
308
309print ("\n" + bcolors.HEADER + "FINISHED!!!" + bcolors.ENDC)