License New BSD license
Lines 413
Keywords
GIS (7) GRASS (1) layer (1) raster (3) tile (1)
Permissions
Viewable by Everyone
Editable by All Siafoo Users
Hide
Solve a problem – Filter by language, license, keyword, owner, or search text to find code & info fast. Join Siafoo Now or Learn More

Layer and Tile GIS Raster Data Atom Feed 0

In Brief Layers and re-tiles GIS rasters using GRASS. Pass in a list of rasters (or directories of rasters), and any null spots in the first raster will be patched in with data from the second, then remaining null spots will be patched with data from the third, and so on. You can specify the desired resolution, extents, and tile size; data will be resampled and sliced appropriately. A color map can also be applied, or raster data can be output as raw. For more information, view the Usage information... more
# 's
  1#!/usr/bin/env python
2
3# Copyright 2007 Regents of the University of California
4# Written by David Isaacson at the University of California, Davis
5# BSD License
6
7import myraster, gdal, getopt, sys, os, gc
8from gdalconst import *
9
10def main():
11 #suppress gdal errors, default handler is 'quiet'
12 gdal.PushErrorHandler()
13
14 paths = []
15 suffixes = ""
16 mapcolor = ""
17 out_command = None
18 output = os.path.abspath(os.getcwd())
19 sizex = sys.maxint
20 sizey = sys.maxint
21 minx = -180
22 miny = -90
23 maxx = 180
24 maxy = 90
25 resolution = 0
26 tiles_only = 0
27 keep_current_tiles = 0
28 force_add = 0
29 delete_inputs = 0
30 delete_outputs = 0
31 raw_output = 0
32
33 my = SystemCommand()
34
35 if not os.getenv("GISBASE"):
36 print "You must be in the GRASS GIS to run this script."
37 sys.exit()
38
39 try:
40 opts,args = getopt.getopt(sys.argv[1:],"a:c:de:ko:qr:s:twx:y:",["command=","force-add","delete-inputs","delete-outputs"])
41 except getopt.GetoptError:
42 #print help information and exit:
43 print "ERROR: Could not parse options."
44 usage()
45 sys.exit(2)
46
47 for opt,arg in opts:
48 if opt == "-a":
49 out_command = arg
50 if opt == "-c":
51 mapcolor = arg
52 if opt == "-d":
53 my.dry_run = 1
54 if opt == "-d":
55 my.dry_run = 1
56 if opt == "-e":
57 extents = arg.split(",")
58 minx, miny, maxx, maxy = \
59 float(extents[0]), float(extents[1]), float(extents[2]), float(extents[3])
60 if opt == "-k":
61 keep_current_tiles = 1
62 if opt == "-o":
63 output = os.path.abspath(arg)
64 if opt == "-q":
65 my.quiet = 1
66 if opt == "-r":
67 resolution = arg
68 if opt == "-s":
69 suffixes = arg.split(",")
70 if opt == "-t":
71 tiles_only = 1
72 if opt == "-w":
73 raw_output = 1
74 if opt == "-x":
75 sizex = long(arg)
76 if opt == "-y":
77 sizey = long(arg)
78 if opt == "--command":
79 args.insert(0,"--command="+arg)
80 if opt == "--force-add":
81 force_add = 1
82 if opt == "--delete-inputs":
83 delete_inputs = 1
84 if opt == "--delete-outputs":
85 delete_outputs = 1
86
87
88 if not len(args):
89 usage()
90 sys.exit(2)
91
92 paths = []
93 command = ""
94 for arg in args:
95 if arg.startswith("--command="):
96 command = arg[10:]
97 command.strip('"')
98 elif os.path.exists(arg):
99 paths.append([arg,command])
100 command = ""
101 else:
102 print arg,"is not a valid path."
103 usage()
104 sys.exit(2)
105
106 if os.path.isdir(output):
107 output = os.path.join(output,"tile")
108
109 if not keep_current_tiles:
110 #find tile sizes
111 print
112 print "Determining tile sizes..."
113 tiles = getTiles(sizex,sizey,minx,miny,maxx,maxy)
114
115 if not tiles:
116 print "ERROR: No valid tiles were found."
117 sys.exit()
118
119 if tiles_only:
120 for tile in tiles:
121 print "TILE (%f,%f) by (%f,%f)" %tuple(tile)
122 sys.exit()
123
124 #make database of extents
125 print
126 print "Looking for rasters and vectors..."
127 dicts = []
128 for rawpath,command in paths:
129 files = myraster.getFilesBySuffix(rawpath,suffixes)
130 for file in files:
131 coords = myraster.findGDALCoordinates(file)
132 is_vector = 0
133 if not coords:
134 coords = myraster.findOGRCoordinates(file)
135 is_vector = 1
136 if not coords: #still
137 continue
138 if not myraster.extentsOverlap([minx,miny,maxx,maxy],coords):
139 continue
140 dict = {'file':file, 'coords':coords, 'raster':os.path.basename(file)}
141 if command:
142 dict['command'] = command
143 if is_vector:
144 vectorName = ""
145 for char in dict['raster']:
146 if char is not ".":
147 vectorName = vectorName + char
148 dict['vector'] = vectorName
149 dicts.append(dict)
150
151 if not dicts:
152 print "ERROR: No valid raster or vector files were found."
153 sys.exit()
154
155 #import rasters
156 if force_add: print "Assuming no rasters or vectors currently exist in grass."
157 print "Importing rasters..."
158
159 existing_rasters = my.system("g.list rast",2)[1].split()
160 existing_vectors = my.system("g.list vect",2)[1].split()
161
162 for dict in dicts:
163 print
164 #raster already exists
165 if dict['raster'] in existing_rasters and not force_add:
166 #if not my.system("g.list rast | grep '%s'" %(dict['raster']),2,[1])[0] and not force_add:
167 print "The raster",dict['raster'],"already exists in grass. Skipping."
168
169 #raster does not exist, import raster
170 elif not dict.has_key('vector'):
171 print "The raster",dict['raster'],"does NOT exist in grass. Importing..."
172 print "Setting region:",
173 my.system("g.region w=%f s=%f e=%f n=%f" %(tuple(dict['coords'])))
174 my.system("r.in.gdal input=%s output=%s" %(dict['file'],dict['raster']))
175
176 #raster does not exist, create from vector
177 else:
178 print "The raster",dict['raster'],"does NOT exist in grass. Creating from vector..."
179 print "Setting region:",
180 my.system("g.region w=%f s=%f e=%f n=%f" %(tuple(dict['coords'])))
181 #work around bug in v.in.ogr --overwrite
182 if force_add: my.system("g.remove vect=%s" %(dict['vector']))
183 #does vector exist?
184 if dict['vector'] in existing_vectors:
185 #if not my.system("g.list vect | grep '%s'" %(dict['vector']),2,[1])[0] and not force_add:
186 print "The vector",dict['vector'],"already exists in grass. Skipping."
187 else:
188 print "The vector",dict['vector'],"does NOT exist in grass. Importing..."
189 vimport = my.system("v.in.ogr -o -e dsn=%s output=%s" %(dict['file'],dict['vector']),handled=[139])[0]
190 if vimport == 139: #something bad happened... we get seg faults occasionally for no apparent reason.
191 print "Vector import failed. Trying alternative approach..."
192 tempvect = dict['vector']+'-temp'
193 my.system("g.remove vect=%s" %(dict['vector']))
194 my.system("v.in.ogr -o -e -c dsn=%s output=%s" %(dict['file'],tempvect))
195 my.system("v.clean input=%s output=%s tool=rmdupl,bpol" %(tempvect,dict['vector']))
196 my.system("g.remove vect=%s" %(tempvect))
197
198 #either way, we want to make our vector into a raster
199 print
200 print "Converting vector to raster..."
201 my.system("v.to.rast --overwrite input=%s output=%s use=val value=-9999" %(dict['vector'],dict['raster']))
202 #now we have our raster, no matter what
203
204 #run command if needed
205 if dict.has_key('command'):
206 command = dict['command']
207 print
208 print "Performing user command on input file..."
209 if command.count("%s") == 1:
210 my.system(command %(dict['raster']))
211 if command.count("%s") == 2:
212 tempname = dict['raster']+'-temp'
213 my.system(command %(dict['raster'],tempname))
214 my.system("g.remove rast=%s" %(dict['raster']))
215 my.system("g.rename rast=%s,%s" %(tempname,dict['raster']))
216 else:
217 my.system(command)
218
219 if not keep_current_tiles:
220 for tile in tiles:
221 print
222 print "******************************"
223 print
224 print "TILE (%f,%f) by (%f,%f)" %tuple(tile)
225
226 rasters = []
227
228 #find rasters for this region
229 print
230 print "Looking for rasters for this tile..."
231 for dict in dicts:
232 raster = dict['raster']
233 coords = dict['coords']
234 if myraster.extentsOverlap(tile,coords):
235 rasters.append(raster)
236 if not rasters:
237 print "No raster files were found for this region. Continuing..."
238 continue
239
240 #set region
241 print
242 print "Setting region..."
243 my.system("g.region w=%f s=%f e=%f n=%f" %tuple(tile))
244
245 #set resolution
246 if resolution:
247 print
248 print "Setting resolution..."
249 my.system("g.region res=%s" %(resolution))
250
251 #get output name
252 outputfile = getFileName(output,tile)
253
254 #create composite
255 print
256 print "Creating composite tile..."
257 composite = os.path.basename(outputfile)
258 if len(rasters) > 1:
259 rasterstring = ",".join(rasters)
260 my.system("r.patch input=%s output=%s --overwrite" %(rasterstring,composite))
261 else:
262 my.system("r.resample input=%s output=%s --overwrite" %(raster,composite))
263
264 #set colors
265 if mapcolor and not raw_output:
266 print
267 print "Setting colors..."
268 my.system("r.colors map=%s color=rules < %s" %(composite,mapcolor))
269
270 #perform output command
271 if out_command:
272 print
273 print "Performing user command on output file..."
274 if out_command.find("%s") != -1:
275 my.system(out_command %(composite))
276 else:
277 my.system(out_command)
278
279 #save to file
280 print
281 print "Saving tile to file..."
282 if raw_output:
283 my.system("r.out.gdal input=%s format=GTiff type=Int16 createopt=\"TILED=YES\" createopt=\"BLOCKXSIZE=128\" createopt=\"BLOCKYSIZE=128\" output=%s" %(composite,outputfile))
284 else:
285 tempbase = output+"-temp"
286 my.system('r.out.tiff -t input="%s" output="%s"' %(composite,tempbase))
287 my.system('gdal_translate -co "TILED=YES" -co "BLOCKXSIZE=128" -co "BLOCKYSIZE=128" -co "COMPRESS=DEFLATE" %s %s ' %(tempbase+".tif",outputfile))
288 my.system('rm -f %s %s' %(tempbase+".tif",tempbase+".tfw"))
289
290 #delete composite
291 print
292 if delete_outputs:
293 print "Deleting composite raster from GRASS..."
294 my.system("g.remove rast=%s" %(composite))
295
296
297 else: #if keep_current_tiles:
298 for dict in dicts:
299 raster = dict['raster']
300 print
301 print "******************************"
302 print
303 print "RASTER:",raster
304
305 #set colors
306 if mapcolor:
307 print
308 print "Setting colors..."
309 #my.system("r.colors rules=%s map=%s" %(mapcolor,raster))
310 my.system("r.colors map=%s color=rules < %s" %(composite,mapcolor))
311
312 if command:
313 print
314 print "Performing given command..."
315 my.system(command %(raster))
316
317 #get output name
318 outputfile = getFileName(output,dict['coords'])
319
320 #save to file
321 print
322 print "Saving raster to file..."
323
324 if raw_output:
325 my.system("r.out.gdal input=%s format=GTiff type=Int16 createopt=\"TILED=YES\" createopt=\"BLOCKXSIZE=128\" createopt=\"BLOCKYSIZE=128\" output=%s" %(raster,outputfile))
326
327 else:
328 tempbase = output+"-temp"
329 my.system("r.out.tiff -t input=%s output=%s" %(raster,tempbase))
330 my.system('gdal_translate -co "TILED=YES,BLOCKXSIZE=128,BLOCKYSIZE=128" %s %s ' %(tempbase+".tif",outputfile))
331 my.system('rm -f %s %s' %(tempbase+".tif",tempbase+".tfw"))
332
333
334 #delete rasters
335 if delete_inputs:
336 print
337 print "Deleting input rasters from GRASS..."
338 for dict in dicts:
339 my.system("g.remove rast=%s" %(dict['raster']))
340
341 print
342 print "myrastertiles.py completed. Check above for errors."
343 print
344
345
346def getTiles(sizex,sizey,minx,miny,maxx,maxy):
347 tiles = []
348 xstripe = getStripes(sizex,minx,maxx)
349 ystripe = getStripes(sizey,miny,maxy)
350 for x in xstripe:
351 for y in ystripe:
352 tiles.append([x[0],y[0],x[1],y[1]])
353 return tiles
354
355
356def getStripes(size,min,max):
357 stripes = []
358 while min+size <= max:
359 stripes.append([min,min+size])
360 min += size
361 if min < max:
362 stripes.append([min,max])
363 return stripes
364
365
366def getFileName(base,coords):
367 if coords[0] >= 0:
368 xstring = "E%03d" %coords[0]
369 else:
370 xstring = "W%03d" %(-1*coords[0])
371 if coords[1] >= 0:
372 ystring = "N%02d" %coords[1]
373 else:
374 ystring = "S%02d" %(-1*coords[1])
375 return base + ystring + xstring + ".tif"
376
377
378class SystemCommand:
379 dry_run = 0
380 quiet = 0
381
382 def system(self,command,quiet=None,handled=[]):
383 #quiet = 0: report command and output
384 #quiet = 1: report command only
385 #quiet = 2: report nothing
386 if quiet is None: quiet=self.quiet
387 if quiet != 2:
388 print command
389 if self.dry_run:
390 return [0,""]
391 command += " 2>&1"
392 fd = os.popen(command,"r",1)
393
394 output = ""
395 while True:
396 char = fd.read(1)
397 if not char:
398 break
399 output += char
400 if quiet == 0:
401 sys.stdout.write(char)
402
403 err = fd.close()
404 if err:
405 err = int(err)/256
406 for handle in handled:
407 if err == handle:
408 return [err,output]
409
410 #on unhandled err, disable quiet
411 if quiet == 2: print command
412 if quiet != 0: sys.stdout.write(output)
413 print
414 print "ERROR: A system call did not exit cleanly."
415 print " Command:",command
416 print " Error Code:",err
417 print "Exiting."
418 print
419 sys.exit()
420
421 return [0,output]
422
423
424def usage():
425 print "Usage:"
426 print " myrastertiles.py [-d] [-k] [-q] [-t] [-w] [-a command] [-c grass_color_file]"
427 print " [-e \"minx,miny,maxx,maxy\"] [-r resolution] [-s suffix1,suffix2,... ]"
428 print " [-x horizontal_size] [-y vertical_size] [-o output_path]"
429 print " [--force-add] [--delete-inputs] [--delete-outputs]"
430 print " [--command=\"g.do.something\"] high_priority_input [ ... "
431 print " [--command=\"g.do.something.else\"] low_priority_input]"
432 print " Inputs may be files or directories."
433 print " Specifying a size smaller than the extent given may result in multiple tiles."
434 print
435 print " Non-obvious options:"
436 print " -a command to run on all output rasters."
437 print " One \"%s\" will be interpreted as the location that the raster name goes."
438 print " Two will be interpreted as input and output raster names, and the output will be"
439 print " copied back over the input."
440 print " Example: \"r.colors map=%s rules=srtm\""
441 print " -d dry run, only print commands to be used."
442 print " -k keep current tiles (-t, -r, -x, -y ignored)"
443 print " -q quiet: do not show the results of system calls"
444 print " -t determine tile widths and resolutions and exit."
445 print " -w output raW (Int16) geoTiff instead of colored (Byte) geoTiff"
446 print " --command command to run on the next input WHEN IMPORTED, same behavior as \"-a\""
447 print " --force-add Always add rasters to GRASS, no matter if"
448 print " appropriately-named rasters already exist."
449 print " --delete-inputs Delete imported rasters from GRASS on exit."
450 print " --delete-outputs Delete created rasters from grass after exporting."
451
452if __name__ == "__main__":
453 main()

Layers and re-tiles GIS rasters using GRASS. Pass in a list of rasters (or directories of rasters), and any null spots in the first raster will be patched in with data from the second, then remaining null spots will be patched with data from the third, and so on. You can specify the desired resolution, extents, and tile size; data will be resampled and sliced appropriately. A color map can also be applied, or raster data can be output as raw. For more information, view the Usage information at the end of the snippet.

Find GIS Files that Overlap Extents / Find Extents of GIS File must be imported (here, as myraster). Incorporates Better Command-Line Calls.

You could probably do most of this in GDAL, but the Python bindings for that are really badly documented, and it would take a while to figure out :-p