[GRASS-SVN] r54581 - in grass-addons/grass7/raster/r.modis: libmodis r.modis.import
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 8 08:03:07 PST 2013
Author: lucadelu
Date: 2013-01-08 08:03:07 -0800 (Tue, 08 Jan 2013)
New Revision: 54581
Modified:
grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py
Log:
fix and cleanup r.modis.import
Modified: grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py 2013-01-08 15:57:56 UTC (rev 54580)
+++ grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py 2013-01-08 16:03:07 UTC (rev 54581)
@@ -64,8 +64,8 @@
# color for lst product
lst_color = ['celsius']
### values of vi product:
- vi_spec = '( 1 1 0 0 0 0 0 0 0 0 0 )'
- vi_specqa = '( 1 1 1 0 0 0 0 0 0 0 0 )'
+ vi_spec = '( 1 1 0 0 0 0 0 0 0 0 0 0 )'
+ vi_specqa = '( 1 1 1 0 0 0 0 0 0 0 0 1 )'
vi_patt = {3: [2, 3], 63: [13, 14, 15], 128: [3], 1024: [1],
8192: [0, 6, 7], 16384: [1], 32768: [1]}
vi_color = ['ndvi', 'evi']
Modified: grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py
===================================================================
--- grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py 2013-01-08 15:57:56 UTC (rev 54580)
+++ grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py 2013-01-08 16:03:07 UTC (rev 54581)
@@ -5,8 +5,8 @@
#
# MODULE: r.modis.import
# AUTHOR(S): Luca Delucchi
-# PURPOSE: r.modis.import is an interface to pyModis for import into GRASS
-# GIS level 3 MODIS produts
+# PURPOSE: r.modis.import is an interface to pyModis for import into
+# GRASS GIS level 3 MODIS produts
#
# COPYRIGHT: (C) 2011 by Luca Delucchi
#
@@ -17,7 +17,7 @@
#############################################################################
#
# REQUIREMENTS:
-# - MRT tools, https://lpdaac.usgs.gov/lpdaac/tools/modis_reprojection_tool
+# - MRT tools, https://lpdaac.usgs.gov/lpdaac/tools/modis_reprojection_tool
#
#############################################################################
@@ -62,7 +62,7 @@
#% key: files
#% type: string
#% key_desc: file
-#% description: Full path to file with list of HDF files
+#% description: Full path to file with list of HDF files
#% gisprompt: old,file,input
#% required: no
#%end
@@ -83,11 +83,15 @@
#% key: spectral
#% type: string
#% key_desc: spectral subset
-#% description: String to choose subset of layers in HDF file for import
+#% description: String to choose subset of layers in HDF file for import
#% required: no
#%end
-import os, sys, string, glob, shutil
+import os
+import sys
+import string
+import glob
+import shutil
import grass.script as grass
from datetime import date
@@ -95,17 +99,18 @@
if os.path.isdir(os.path.join(os.getenv('GISBASE'), 'etc', 'r.modis')):
libmodis = os.path.join(os.getenv('GISBASE'), 'etc', 'r.modis')
elif os.getenv('GRASS_ADDON_BASE') and \
- os.path.isdir(os.path.join(os.getenv('GRASS_ADDON_BASE'), 'etc', 'r.modis')):
+ os.path.isdir(os.path.join(os.getenv('GRASS_ADDON_BASE'), 'etc',
+ 'r.modis')):
libmodis = os.path.join(os.getenv('GRASS_ADDON_BASE'), 'etc', 'r.modis')
elif os.path.isdir(os.path.join('..', 'libmodis')):
- libmodis = os.path.join('..', 'libmodis')
+ libmodis = os.path.join('..', 'libmodis')
if not libmodis:
sys.exit("ERROR: modis library not found")
sys.path.append(libmodis)
# try to import pymodis (modis) and some classes for r.modis.download
try:
- from rmodislib import resampling, product, get_proj, projection
+ from rmodislib import resampling, product, projection
except ImportError, e:
grass.fatal(e)
try:
@@ -117,13 +122,15 @@
except ImportError, e:
grass.fatal(e)
-def list_files(opt, mosaik = False):
- """If used in function single(): Return a list of HDF files from the file list
- If used in function mosaic(): Return a dictionary with a list of HDF files for each day
+
+def list_files(opt, mosaik=False):
+ """If used in function single(): Return a list of HDF files from the file
+ list. If used in function mosaic(): Return a dictionary with a list of HDF
+ files for each day
"""
# read the file with the list of HDF
if opt['files'] != '':
- listoffile = open(opt['files'],'r')
+ listoffile = open(opt['files'], 'r')
basedir = os.path.split(listoffile.name)[0]
# if mosaic create a dictionary
if mosaik:
@@ -133,10 +140,10 @@
filelist = []
# append hdf files
for line in listoffile:
- if string.find(line,'xml') == -1 and mosaik == False:
+ if string.find(line, 'xml') == -1 and mosaik == False:
filelist.append(line.strip())
# for mosaic create a list of hdf files for each day
- elif string.find(line,'xml') == -1 and mosaik == True:
+ elif string.find(line, 'xml') == -1 and mosaik == True:
day = line.split('/')[-1].split('.')[1]
if filelist.has_key(day):
filelist[day].append(line)
@@ -148,7 +155,8 @@
basedir = os.path.split(filelist[0])[0]
return filelist, basedir
-def spectral(opts, prod, q, m = False):
+
+def spectral(opts, prod, q, m=False):
"""Return spectral string"""
# return the spectral set selected by the user
if opts['spectral'] != '':
@@ -158,7 +166,7 @@
if q:
if prod['spec_qa']:
spectr = prod['spec_qa']
- else:
+ else:
spectr = prod['spec']
else:
spectr = prod['spec']
@@ -166,9 +174,10 @@
spectr = spectr.replace(' 0', '')
return spectr
+
def confile(pm, opts, q, mosaik=False):
"""Create the configuration file for MRT software"""
- # return projection and datum
+ # return projection and datum
projObj = projection()
proj = projObj.returned()
dat = projObj.datum()
@@ -179,7 +188,7 @@
cod = os.path.split(pm.hdfname)[1].split('.')[0]
prod = product().fromcode(cod)
if mosaik:
- # if mosaic it remove all the 0 from the subset string to convert all
+ # if mosaic it remove all the 0 from the subset string to convert all
# the right layer
spectr = spectral(opts, prod, q, True)
else:
@@ -198,26 +207,30 @@
else:
res = None
try:
- conf = pm.confResample(spectr, res, pref, dat, resampl,
+ conf = pm.confResample(spectr, res, pref, dat, resampl,
proj, zone, projpar)
return conf
except IOError, e:
grass.fatal(e)
-def prefix(options, name = False):
+
+def prefix(options, name=False):
"""Return the prefix of output file if not set return None to use default
value
"""
if options['output']:
- return options['output']
+ return '%s.tif' % options['output']
else:
- return None
+ return None
+
def import_tif(out, basedir, rem, write, target=None):
"""Import TIF files"""
# list of tif files
- tifiles = glob.glob1(basedir, out + "*.tif")
+ tifiles = glob.glob1(basedir, "*.tif")
if not tifiles:
+ tifiles = glob.glob1(os.getcwd(), "*.tif")
+ if not tifiles:
grass.fatal(_('Error during the conversion'))
# check if user is in latlong location to set flag l
if projection().val == 'll':
@@ -228,25 +241,32 @@
# for each file import it
for t in tifiles:
basename = os.path.splitext(t)[0]
- name = os.path.join(basedir,t)
+ name = os.path.join(basedir, t)
+ if not os.path.exists(name):
+ name = os.path.join(os.getcwd(), t)
+ if not os.path.exists(name):
+ grass.warning(_("File %s doesn't find" % basename))
+ continue
filesize = int(os.path.getsize(name))
if filesize < 1000:
grass.warning(_('Probably some error occur during the conversion' \
+ 'for file <%s>. Escape import' % name))
continue
try:
- grass.run_command('r.in.gdal', input = name, output = basename,
- flags = f, overwrite = write, quiet = True)
+ grass.run_command('r.in.gdal', input=name, output=basename,
+ flags=f, overwrite=write, quiet=True)
outfile.append(basename)
except:
- grass.fatal(_('Error during import'))
+ grass.warning(_('Error during import of %s' % basename))
+ continue
if rem:
os.remove(name)
- if target:
+ if target:
if target != basedir:
- shutil.move(name,target)
+ shutil.move(name, target)
return outfile
+
def findfile(pref, suff):
""" Check if a file exists on mapset """
if grass.find_file(pref + suff)['file']:
@@ -254,32 +274,38 @@
else:
grass.warning(_("Raster map <%s> not found") % (pref + suff))
+
def metadata(pars, mapp, coll):
""" Set metadata to the imported files """
# metadata
meta = pars.metastring()
- grass.run_command('r.support', quiet = True, map=mapp, hist=meta)
+ grass.run_command('r.support', quiet=True, map=mapp, hist=meta)
# timestamp
rangetime = pars.retRangeTime()
data = rangetime['RangeBeginningDate'].split('-')
- dataobj = date(int(data[0]),int(data[1]),int(data[2]))
- grass.run_command('r.timestamp', map=mapp, date=dataobj.strftime("%d %b %Y"),
- quiet = True)
+ dataobj = date(int(data[0]), int(data[1]), int(data[2]))
+ grass.run_command('r.timestamp', map=mapp, quiet=True,
+ date=dataobj.strftime("%d %b %Y"))
# color
if string.find(mapp, 'QC') != -1 or string.find(mapp, 'Quality') != -1 or \
string.find(mapp, 'QA') != -1:
- grass.run_command('r.colors', quiet = True, map=mapp, color = coll)
+ grass.run_command('r.colors', quiet=True, map=mapp, color=coll)
elif string.find(mapp, 'NDVI') != -1:
- grass.run_command('r.colors', quiet = True, map=mapp, color = coll[0])
+ grass.run_command('r.colors', quiet=True, map=mapp, color=coll[0])
elif string.find(mapp, 'EVI') != -1:
- grass.run_command('r.colors', quiet = True, map=mapp, color = coll[1])
+ grass.run_command('r.colors', quiet=True, map=mapp, color=coll[1])
elif string.find(mapp, 'LST') != -1:
- grass.run_command('r.colors', quiet = True, map=mapp, color = coll[0])
+ grass.run_command('r.colors', quiet=True, map=mapp, color=coll[0])
elif string.find(mapp, 'Snow') != -1:
- grass.run_command('r.colors', quiet = True, map=mapp, color = coll[0])
+ grass.run_command('r.colors', quiet=True, map=mapp, color=coll[0])
+
def analyze(pref, an, cod, parse, write):
""" Analyze the MODIS data using QA if present """
+ if pref.find('.tif') != -1:
+ pref = pref.rstrip('.tif')
+ import pdb
+ pdb.set_trace()
prod = product().fromcode(cod)
if not prod['spec_qa']:
grass.warning(_("There is no QA layer, analysis and filtering will be skipped"))
@@ -289,113 +315,124 @@
col = prod['color']
val = []
qa = []
- for v,q in suf.iteritems():
- val.append(findfile(pref,v))
+ for v, q in suf.iteritems():
+ val.append(findfile(pref, v))
if q:
- qa.append(findfile(pref,q))
+ qa.append(findfile(pref, q))
+ pdb.set_trace()
for n in range(len(val)):
if val[n] == None:
grass.warning(_("Some error occur"))
continue
valname = val[n]['name']
valfull = val[n]['fullname']
- grass.run_command('g.region', rast = valfull)
- grass.run_command('r.null', map = valfull, setnull = 0)
- if string.find(cod,'13Q1') >= 0 or string.find(cod,'13A2') >= 0:
+ grass.run_command('g.region', rast=valfull)
+ grass.run_command('r.null', map=valfull, setnull=0)
+ if string.find(cod, '13Q1') >= 0 or string.find(cod, '13A2') >= 0:
mapc = "%s.2 = %s / 10000." % (valname, valfull)
- elif string.find(cod,'11A1') >= 0 or string.find(cod,'11A2') >= 0 or string.find(cod,'11B1') >= 0:
+ elif string.find(cod, '11A1') >= 0 or string.find(cod, '11A2') >= 0 \
+ or string.find(cod, '11B1') >= 0:
mapc = "%s.2 = (%s * 0.0200) - 273.15" % (valname, valfull)
+ pdb.set_trace()
grass.mapcalc(mapc)
if an == 'noqa':
- #grass.run_command('g.remove', quiet = True, rast = valfull)
+ #grass.run_command('g.remove', quiet=True, rast = valfull)
try:
- grass.run_command('g.rename', quiet = True, overwrite = write,
- rast=(valname,valname + '.orig'))
- grass.run_command('g.rename', quiet = True, overwrite = write,
- rast= (valname + '.2', valname))
+ grass.run_command('g.rename', quiet=True, overwrite=write,
+ rast=(valname, valname + '.orig'))
+ grass.run_command('g.rename', quiet=True, overwrite=write,
+ rast=(valname + '.2', valname))
except:
pass
metadata(parse, valname, col)
metadata(parse, valname, col)
- metadata(parse, qafull, 'byr')
+ metadata(parse, valname, 'byr')
if an == 'all':
if len(qa) != len(val):
- grass.fatal(_("The number of QA and value maps is different, something is wrong"))
+ grass.fatal(_("The number of QA and value maps is different,"\
+ " something is wrong"))
qaname = qa[n]['name']
qafull = qa[n]['fullname']
finalmap = "%s.3=if(" % valname
first_map = 1
- for key,value in prod['pattern'].iteritems():
+ for key, value in prod['pattern'].iteritems():
for v in value:
outpat = "%s.%i.%i" % (qaname, key, v)
- grass.run_command('r.bitpattern', quiet = True, input = qafull,
- output = outpat, pattern = key, patval= v)
+ grass.run_command('r.bitpattern', quiet=True, input=valname,
+ output=outpat, pattern=key, patval=v)
if first_map:
first_map = 0
finalmap += "%s == 0 " % outpat
else:
finalmap += "&& %s == 0 " % outpat
- if string.find(cod,'13Q1') >= 0 or string.find(cod,'13A2') >= 0:
+ if string.find(cod, '13Q1') >= 0 or string.find(cod, '13A2') >= 0:
finalmap += "&& %s.2 <= 1.000" % valname
finalmap += ",%s.2, null() )" % valname
# grass.message("mapc finalmap: %s" % finalmap)
+ pdb.set_trace()
grass.mapcalc(finalmap)
- #grass.run_command('g.remove', quiet = True, rast=(valname, valname + '.2'))
- grass.run_command('g.rename', quiet = True, overwrite = write,
- rast=(valname,valname + '.orig'))
- grass.run_command('g.remove', quiet = True, rast=(valname + '.2'))
- grass.run_command('g.mremove', flags="f", quiet = True,
- rast = ("%s.*" % qaname))
- grass.run_command('g.rename', quiet = True, overwrite = write,
+ grass.run_command('g.rename', quiet=True, overwrite=write,
+ rast=(valname, valname + '.orig'))
+ grass.run_command('g.remove', quiet=True, rast=(valname + '.2'))
+ grass.run_command('g.mremove', flags="f", quiet=True,
+ rast=("%s.*" % qaname))
+ grass.run_command('g.rename', quiet=True, overwrite=write,
rast=(valname + '.3', valname))
metadata(parse, valname, col)
metadata(parse, valname, col)
- metadata(parse, qafull, 'byr')
+ metadata(parse, valname, 'byr')
-def single(options,remove,an,ow):
+
+def single(options, remove, an, ow):
"""Convert the HDF file to TIF and import it
"""
listfile, basedir = list_files(options)
pid = str(os.getpid())
# for each file
for i in listfile:
- # the full path to hdf file
- hdf = os.path.join(basedir,i)
+ if os.path.exists(i):
+ hdf = i
+ else:
+ # the full path to hdf file
+ hdf = os.path.join(basedir, i)
+ if not os.path.exists(hdf):
+ grass.warning(_("%s not found" % i))
+ continue
# create conf file fro mrt tools
pm = parseModis(hdf)
- confname = confile(pm,options,an)
- print confname
+ confname = confile(pm, options, an)
# create convertModis class and convert it in tif file
- execmodis = convertModis(hdf,confname,options['mrtpath'])
+ execmodis = convertModis(hdf, confname, options['mrtpath'])
execmodis.run()
output = prefix(options)
if not output:
output = os.path.split(hdf)[1].rstrip('.hdf')
# import tif files
- maps_import = import_tif(output,basedir,remove,ow)
+ maps_import = import_tif(output, basedir, remove, ow)
if an and len(maps_import) != 0:
- grass.run_command('g.region', save = 'oldregion.%s' % pid)
+ grass.run_command('g.region', save='oldregion.%s' % pid)
try:
cod = os.path.split(pm.hdfname)[1].split('.')[0]
- analyze(outname, an, cod, pm, ow)
+ analyze(output, an, cod, pm, ow)
except:
- grass.run_command('g.region', region = 'oldregion.%s' % pid)
- grass.run_command('g.remove',quiet = True,
- region = 'oldregion.%s' % pid)
- cod = os.path.split(pm.hdfname)[1].split('.')[0]
- analyze(output, an, cod, pm, ow)
+ grass.run_command('g.region', region='oldregion.%s' % pid)
+ grass.run_command('g.remove', quiet=True,
+ region='oldregion.%s' % pid)
+# cod = os.path.split(pm.hdfname)[1].split('.')[0]
+# analyze(output, an, cod, pm, ow)
#os.remove(confname)
-def mosaic(options,remove,an,ow):
+
+def mosaic(options, remove, an, ow):
"""Create a daily mosaic of HDF files convert to TIF and import it
"""
- dictfile, targetdir = list_files(options,True)
+ dictfile, targetdir = list_files(options, True)
pid = str(os.getpid())
# for each day
for dat, listfiles in dictfile.iteritems():
# create the file with the list of name
- tempfile = open(os.path.join(targetdir,pid),'w')
+ tempfile = open(os.path.join(targetdir, pid), 'w')
tempfile.writelines(listfiles)
tempfile.close()
# basedir of tempfile, where hdf files are write
@@ -404,7 +441,7 @@
listfiles[0].split('/')[-1].split('.')[1])
# return the spectral subset in according mrtmosaic tool format
prod = product().fromcode(listfiles[0].split('/')[-1].split('.')[0])
- spectr = spectral(options,prod,an)
+ spectr = spectral(options, prod, an)
spectr = spectr.lstrip('( ').rstrip(' )')
# create mosaic
cm = createMosaic(tempfile.name, outname, options['mrtpath'], spectr)
@@ -413,26 +450,26 @@
hdfiles = glob.glob1(basedir, outname + "*.hdf")
for i in hdfiles:
# the full path to hdf file
- hdf = os.path.join(basedir,i)
+ hdf = os.path.join(basedir, i)
# create conf file fro mrt tools
pm = parseModis(hdf)
- confname = confile(pm,options,an, True)
+ confname = confile(pm, options, an, True)
# create convertModis class and convert it in tif file
execmodis = convertModis(hdf, confname, options['mrtpath'])
execmodis.run()
- # remove hdf
- if remove:
+ # remove hdf
+ if remove:
# import tif files
maps_import = import_tif(outname, basedir, remove, ow)
if an:
- grass.run_command('g.region', save = 'oldregion.%s' % pid)
+ grass.run_command('g.region', save='oldregion.%s' % pid)
try:
cod = os.path.split(pm.hdfname)[1].split('.')[0]
analyze(outname, an, cod, pm, ow)
except:
- grass.run_command('g.region', region = 'oldregion.%s' % pid)
- grass.run_command('g.remove',quiet = True,
- region = 'oldregion.%s' % pid)
+ grass.run_command('g.region', region='oldregion.%s' % pid)
+ grass.run_command('g.remove', quiet=True,
+ region='oldregion.%s' % pid)
os.remove(hdf)
os.remove(hdf + '.xml')
# or move the hdf and hdf.xml to the dir where are the original files
@@ -440,23 +477,24 @@
# import tif files
import_tif(outname, basedir, remove, ow, targetdir)
if an and len(maps_import) != 0:
- grass.run_command('g.region', save = 'oldregion.%s' % pid)
+ grass.run_command('g.region', save='oldregion.%s' % pid)
try:
cod = os.path.split(pm.hdfname)[1].split('.')[0]
analyze(outname, an, cod, pm, ow)
except:
- grass.run_command('g.region', region = 'oldregion.%s' % pid)
- grass.run_command('g.remove', region = 'oldregion.%s' % pid)
- try:
- shutil.move(hdf,targetdir)
- shutil.move(hdf + '.xml',targetdir)
+ grass.run_command('g.region', region='oldregion.%s' % pid)
+ grass.run_command('g.remove', region='oldregion.%s' % pid)
+ try:
+ shutil.move(hdf, targetdir)
+ shutil.move(hdf + '.xml', targetdir)
except:
pass
# remove the conf file
os.remove(confname)
grass.try_remove(tempfile.name)
- grass.try_remove(os.path.join(targetdir,'mosaic',pid))
+ grass.try_remove(os.path.join(targetdir, 'mosaic', pid))
+
def main():
# check if you are in GRASS
gisbase = os.getenv('GISBASE')
@@ -501,9 +539,9 @@
grass.fatal(_('It is not possible to create a mosaic with a single HDF file'))
return 0
elif flags['m']:
- mosaic(options,remove,analyze,over)
+ mosaic(options, remove, analyze, over)
else:
- single(options,remove,analyze,over)
+ single(options, remove, analyze, over)
if __name__ == "__main__":
options, flags = grass.parser()
More information about the grass-commit
mailing list