3
############################################################################
5
# MODULE: r_in_aster.py
6
# AUTHOR(S): Michael Barton (michael.barton@asu.edu) and
7
# Glynn Clements (glynn@gclements.plus.com)
8
# Based on r.in.aster bash script for GRASS
9
# by Michael Barton and Paul Kelly
10
# PURPOSE: Rectifies, georeferences, & imports Terra-ASTER imagery
12
# COPYRIGHT: (C) 2008 by the GRASS Development Team
14
# This program is free software under the GNU General Public
15
# License (>=v2). Read the file COPYING that comes with GRASS
18
#############################################################################
22
# gdal compiled with HDF4 support
25
#% description: Georeference, rectify, and import Terra-ASTER imagery and relative DEMs using gdalwarp.
29
#% keyword: Terra-ASTER
31
#%option G_OPT_F_INPUT
32
#% description: Name of input ASTER image
37
#% description: ASTER imagery processing type (Level 1A, Level 1B, or relative DEM)
38
#% options: L1A,L1B,DEM
45
#% description: List L1A or L1B band to translate (1,2,3n,...), or enter 'all' to translate all bands
49
#%option G_OPT_R_OUTPUT
50
#% description: Base name for output raster map (band number will be appended to base name)
56
import grass.script as grass
60
'1': "VNIR_Band1:ImageData",
61
'2': "VNIR_Band2:ImageData",
62
'3n': "VNIR_Band3N:ImageData",
63
'3b': "VNIR_Band3B:ImageData",
64
'4': "SWIR_Band4:ImageData",
65
'5': "SWIR_Band5:ImageData",
66
'6': "SWIR_Band6:ImageData",
67
'7': "SWIR_Band7:ImageData",
68
'8': "SWIR_Band8:ImageData",
69
'9': "SWIR_Band9:ImageData",
70
'10': "TIR_Band10:ImageData",
71
'11': "TIR_Band11:ImageData",
72
'12': "TIR_Band12:ImageData",
73
'13': "TIR_Band13:ImageData",
74
'14': "TIR_Band14:ImageData"
77
'1': "VNIR_Swath:ImageData1",
78
'2': "VNIR_Swath:ImageData2",
79
'3n': "VNIR_Swath:ImageData3N",
80
'3b': "VNIR_Swath:ImageData3B",
81
'4': "SWIR_Swath:ImageData4",
82
'5': "SWIR_Swath:ImageData5",
83
'6': "SWIR_Swath:ImageData6",
84
'7': "SWIR_Swath:ImageData7",
85
'8': "SWIR_Swath:ImageData8",
86
'9': "SWIR_Swath:ImageData9",
87
'10': "TIR_Swath:ImageData10",
88
'11': "TIR_Swath:ImageData11",
89
'12': "TIR_Swath:ImageData12",
90
'13': "TIR_Swath:ImageData13",
91
'14': "TIR_Swath:ImageData14"
96
input = options['input']
97
proctype = options['proctype']
98
output = options['output']
99
band = options['band']
101
#check whether gdalwarp is in path and executable
102
if not grass.find_program('gdalwarp', '--help'):
103
grass.fatal(_("gdalwarp is not in the path and executable"))
105
#create temporary file to hold gdalwarp output before importing to GRASS
106
tempfile = grass.read_command("g.tempfile", pid = os.getpid()).strip() + '.tif'
108
#get projection information for current GRASS location
109
proj = grass.read_command('g.proj', flags = 'jf').strip()
111
#currently only runs in projected location
112
if "XY location" in proj:
113
grass.fatal(_("This module needs to be run in a projected location (found: %s)") % proj)
116
#process list of bands
117
allbands = ['1','2','3n','3b','4','5','6','7','8','9','10','11','12','13','14']
121
bandlist = band.split(',')
123
#initialize datasets for L1A and L1B
124
if proctype in ["L1A", "L1B"]:
125
for band in bandlist:
127
dataset = bands[proctype][band]
128
srcfile = "HDF4_EOS:EOS_SWATH:%s:%s" % (input, dataset)
129
import_aster(proj, srcfile, tempfile, band)
131
grass.fatal(_('band %s is not an available Terra/ASTER band') % band)
132
elif proctype == "DEM":
134
import_aster(proj, srcfile, tempfile, "DEM")
137
grass.message(_("Cleaning up ..."))
138
grass.try_remove(tempfile)
139
grass.message(_("Done."))
143
def import_aster(proj, srcfile, tempfile, band):
144
#run gdalwarp with selected options (must be in $PATH)
145
#to translate aster image to geotiff
146
grass.message(_("Georeferencing aster image ..."))
147
grass.debug("gdalwarp -t_srs %s %s %s" % (proj, srcfile, tempfile))
149
if platform.system() == "Darwin":
150
cmd = ["arch", "-i386", "gdalwarp", "-t_srs", proj, srcfile, tempfile ]
152
cmd = ["gdalwarp", "-t_srs", proj, srcfile, tempfile ]
155
#check to see if gdalwarp executed properly
158
#import geotiff to GRASS
159
grass.message(_("Importing into GRASS ..."))
160
outfile = "%s.%s" % (output, band)
161
grass.run_command("r.in.gdal", input = tempfile, output = outfile)
164
grass.raster_history(outfile)
166
if __name__ == "__main__":
167
options, flags = grass.parser()