[GRASS-SVN] r47639 - in grass-addons/grass7/raster/r.modis:
libmodis r.modis.import
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Aug 15 03:56:40 EDT 2011
Author: lucadelu
Date: 2011-08-15 00:56:40 -0700 (Mon, 15 Aug 2011)
New Revision: 47639
Modified:
grass-addons/grass7/raster/r.modis/libmodis/modis.py
grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py
Log:
some fix, add functionality to analize modis data with qa map
Modified: grass-addons/grass7/raster/r.modis/libmodis/modis.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/modis.py 2011-08-14 17:55:15 UTC (rev 47638)
+++ grass-addons/grass7/raster/r.modis/libmodis/modis.py 2011-08-15 07:56:40 UTC (rev 47639)
@@ -49,6 +49,7 @@
import socket
from ftplib import FTP
import ftplib
+from xml.etree import ElementTree
class downModis:
"""A class to download MODIS data from NASA FTP repository"""
@@ -406,7 +407,6 @@
"""Initialization function :
filename = the name of MODIS hdf file
"""
- from xml.etree import ElementTree
if os.path.exists(filename):
# hdf name
self.hdfname = filename
@@ -530,8 +530,9 @@
self.getGranule()
rangeTime = {}
for i in self.granule.find('RangeDateTime').getiterator():
- if i.text.strip() != '':
- rangeTime[i.tag] = i.text
+ if i.text != None:
+ if i.text.strip() != '':
+ rangeTime[i.tag] = i.text
return rangeTime
def retBoundary(self):
@@ -611,6 +612,19 @@
self.getGranule()
return self.granule.find('BrowseProduct').find('BrowseGranuleId').text
+ def metastring(self):
+ date = self.retRangeTime()
+ qa = self.retMeasure()
+ pge = self.retPGEVersion()
+ daynight = self.retDataGranule()
+ stri = "Date:%s;BeginningTime:%s;EndingTime:%s;DayNightFlag:%s;"
+ stri += "QAPercentCloudCover:%s;QAPercentMissingData:%s;PGEVersion:%s"
+ out = stri % ( date['RangeBeginningDate'], date['RangeBeginningTime'],
+ date['RangeEndingTime'], daynight['DayNightFlag'],
+ qa['QAStats']['QAPercentCloudCover'],
+ qa['QAStats']['QAPercentMissingData'], pge)
+ return out
+
def confResample(self, spectral, res = None, output = None, datum = 'WGS84',
resampl = 'NEAREST_NEIGHBOR', projtype = 'GEO', utm = None,
projpar = '( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )',
@@ -786,7 +800,6 @@
"""A class to some variable for the xml file of a mosaic
"""
def __init__(self,hdflist):
- from xml.etree import ElementTree
self.ElementTree = ElementTree
self.hdflist = hdflist
self.parModis = []
@@ -839,38 +852,14 @@
elem = self.ElementTree.SubElement(ele,k)
elem.text = v
- def valDTD(self,obj):
+ def moreValues(self, fun, obj, ele):
values = []
for i in self.parModis:
- values.append(i.retDTD())
+ values.append(getattr(i,fun)())
for i in self._checkval(values):
- dtd = self.ElementTree.SubElement(obj,'DTDVersion')
+ dtd = self.ElementTree.SubElement(obj,ele)
dtd.text = i
- def valDataCenter(self,obj):
- values = []
- for i in self.parModis:
- values.append(i.retDataCenter())
- for i in self._checkval(values):
- dci = self.ElementTree.SubElement(obj,'DataCenterId')
- dci.text = i
-
- def valGranuleUR(self,obj):
- values = []
- for i in self.parModis:
- values.append(i.retGranuleUR())
- for i in self._checkval(values):
- gur = self.ElementTree.SubElement(obj,'GranuleUR')
- gur.text = i
-
- def valDbID(self,obj):
- values = []
- for i in self.parModis:
- values.append(i.retDbID())
- for i in self._checkval(values):
- dbid = self.ElementTree.SubElement(obj,'DbID')
- dbid.text = i
-
def valInsTime(self,obj):
values = []
for i in self.parModis:
@@ -933,7 +922,6 @@
def valInputPointer(self,obj):
for i in self.parModis:
- import pdb; pdb.set_trace()
for v in i.retInputGranule():
ip = self.ElementTree.SubElement(obj,'InputPointer')
ip.text = v
@@ -976,16 +964,17 @@
# the root element
granule = self.ElementTree.Element('GranuleMetaDataFile')
# add DTDVersion
- self.valDTD(granule)
+ self.moreValues('retDTD',granule,'DTDVersion')
+ #self.valDTD(granule)
# add DataCenterId
- self.valDataCenter(granule)
+ self.moreValues('retDataCenter',granule,'DataCenterId')
# add GranuleURMetaData
gurmd = self.ElementTree.SubElement(granule,'GranuleURMetaData')
# add GranuleUR
- self.valGranuleUR(gurmd)
+ self.moreValues('retGranuleUR',granule,'GranuleUR')
# add dbID
- self.valDbID(gurmd)
-
+ self.moreValues('retDbID',granule,'DbID')
+
# TODO ADD InsertTime LastUpdate
# add CollectionMetaData
@@ -1129,17 +1118,20 @@
self.subset = subset
def write_mosaic_xml(self):
+ self.finalfile = open(os.path.join(self.basepath,'mosaic%i' % os.getpid()),'w')
listHDF = []
for i in self.HDFfiles:
if i.find(self.basepath) == -1:
- print "Attection maybe you have the not full path in the HDF file list"
+ print "Attection maybe you are not using full path in the HDF file list"
listHDF.append(os.path.join(self.basepath,i.strip()))
+ self.finalfile.write("%s\n" % os.path.join(self.basepath,i.strip()))
else:
listHDF.append(i.strip())
+ self.finalfile.write("%s\n" % i.strip())
+ self.finalfile.close()
pmm = parseModisMulti(listHDF)
pmm.writexml(self.outxml)
-
def executable(self):
"""Return the executable of mrtmosaic MRT software
"""
@@ -1160,11 +1152,12 @@
else:
self.write_mosaic_xml()
if self.subset:
- subprocess.call([execut,'-i',self.listfiles,'-o',self.out,'-s',self.subset],
+ subprocess.call([execut,'-i',self.finalfile.name,'-o',self.out,'-s',self.subset],
stderr = subprocess.STDOUT)
else:
- subprocess.call([execut,'-i',self.listfiles,'-o',self.out], stderr =
+ subprocess.call([execut,'-i',self.finalfile.name,'-o',self.out], stderr =
subprocess.STDOUT)
+ #os.remove(self.finalfile)
return "The mosaic file %s is created" % self.out
class processModis:
@@ -1225,4 +1218,4 @@
+ 'software') % execut
else:
subprocess.call([execut,'-pf=%s' % self.conf])
- return "The hdf file %s was converted" % self.name
\ No newline at end of file
+ return "The hdf file %s was converted" % self.name
Modified: grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py 2011-08-14 17:55:15 UTC (rev 47638)
+++ grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py 2011-08-15 07:56:40 UTC (rev 47639)
@@ -15,7 +15,7 @@
# for details.
#
#############################################################################
-from grass.script import read_command,parse_command
+import grass.script as grass
# interface to g.proj -p
def get_proj():
"""!Returns the output from running "g.proj -p" plus towgs84 parameter (g.proj -d),
@@ -29,7 +29,7 @@
@return dictionary of projection values
"""
- gproj = read_command('g.proj',flags='p')
+ gproj = grass.read_command('g.proj',flags='p')
listproj = gproj.split('\n')
listproj.remove('-PROJ_INFO-------------------------------------------------')
listproj.remove('-PROJ_UNITS------------------------------------------------')
@@ -38,7 +38,7 @@
for i in listproj:
ilist = i.split(':')
proj[ilist[0].strip()] = ilist[1].strip()
- proj.update(parse_command('g.proj',flags='j'))
+ proj.update(grass.parse_command('g.proj',flags='j'))
return proj
class product:
@@ -47,22 +47,93 @@
def __init__(self,value = None):
urlbase = 'e4ftl01u.ecs.nasa.gov'
usrsnow = 'n4ftl01u.ecs.nasa.gov'
+ lst_spec = '( 1 0 0 0 1 0 0 0 0 0 0 0 )'
+ lst_specqa = '( 1 1 0 0 1 1 0 0 0 0 0 0 )'
+ lstvi_specall = '( 1 1 1 1 1 1 1 1 1 1 1 1)'
+ lst_patt = [ " == 2", " == 3", " >= 81" ]
+ lst_suff = {'.LST_Day_1km':'.QC_Day','.LST_Night_1km':'.QC_Night'}
+ lst_color = ['modis_lst']
+ 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_patt = []
+ vi_color = ['ndvi','evi']
+ vi_suff = {'.250m_16_days_NDVI.tif' : '250m_16_days_VI_Quality.tif',
+ '.250m_16_days_EVI.tif' : '250m_16_days_VI_Quality.tif'}
+
+ snow_spec = ('( 1 1 )')
+ snow_color = ('evi') #TODO CREATE THE COLOR FOR MODIS_SNOW
+ lstL2_spec = 'LST; QC; Error_LST; Emis_31; Emis_32; View_angle; View_time'
+
self.prod = value
- self.products = {
- 'lst_aqua_daily':{'url':urlbase,'folder':'MOLA/MYD11A1.005','res':1000,
- 'spec':'( 1 0 0 0 1 0 0 0 0 0 0 0 )','spec_qa':'( 1 1 0 0 1 1 0 0 0 0 0 0 )'},
- 'lst_terra_daily':{'url':urlbase,'folder':'MOLT/MOD11A1.005','res':1000,
- 'spec':'( 1 0 0 0 1 0 0 0 0 0 0 0 )','spec_qa':'( 1 1 0 0 1 1 0 0 0 0 0 0 )'},
- 'snow_terra_eight':{'url':usrsnow,'folder':'SAN/MOST/MOD10A2.005','res':500,
- 'spec':'( 1 1 0 0 0 0 0 0 0 0 0 )','spec_qa':'( 0 0 0 0 0 0 0 0 0 0 0 )'},
- 'ndvi_terra_sixte':{'url':urlbase, 'folder':'MOLT/MOD13Q1.005','res':250,
- 'spec':'( 1 1 0 0 0 0 0 0 0 0 0 )','spec_qa':'( 1 1 1 1 0 0 0 0 0 0 0 )'}
+ lst = {'lst_aqua_daily_1000' : {'url' : urlbase, 'folder' : 'MOLA/MYD11A1.005',
+ 'res' : 1000, 'spec' : lst_spec, 'spec_qa' : lst_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : lst_suff,
+ 'pattern' : lst_patt, 'color' : lst_color},
+ 'lst_terra_daily_1000' : {'url' : urlbase, 'folder' : 'MOLT/MOD11A1.005',
+ 'res' : 1000, 'spec': lst_spec,'spec_qa': lst_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : lst_suff,
+ 'pattern' : lst_patt, 'color' : lst_color},
+ 'lst_terra_eight_1000' : {'url' : urlbase, 'folder' : 'MOLT/MOD11A2.005',
+ 'res' : 1000, 'spec': lst_spec,'spec_qa': lst_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : lst_suff,
+ 'pattern' : lst_patt, 'color' : lst_color},
+ 'lst_aqua_eight_1000' : {'url' : urlbase, 'folder' : 'MOLA/MYD11A2.005',
+ 'res' : 1000, 'spec': lst_spec,'spec_qa': lst_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : lst_suff,
+ 'pattern' : lst_patt, 'color' : lst_color},
+ 'lst_terra_daily_6000' : {'url' : urlbase, 'folder' : 'MOLT/MOD11B1.005',
+ 'res' : 6000, 'spec': lst_spec,'spec_qa': lst_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : lst_suff,
+ 'pattern' : lst_patt, 'color' : lst_color},
+ 'lst_aqua_daily_6000' : {'url' : urlbase, 'folder' : 'MOLA/MYD11B1.005',
+ 'res' : 6000, 'spec': lst_spec,'spec_qa': lst_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : lst_suff,
+ 'pattern' : lst_patt, 'color' : lst_color},
+
}
- self.products_swath = {
- 'lstemi_terra_daily':{'url':urlbase,'folder':'MOLT/MOD11_L2.005',
- 'spec':'LST; QC; Error_LST; Emis_31; Emis_32; View_angle; View_time'}, 'lstemi_aqua_daily':{'url':urlbase,'folder':'MOLA/MYD11_L2.005',
- 'spec':'LST; QC; Error_LST; Emis_31; Emis_32; View_angle; View_time'}
+ vi = {'ndvi_terra_sixteen_250':{'url':urlbase, 'folder':'MOLT/MOD13Q1.005',
+ 'res':250,'spec': vi_spec,'spec_qa': vi_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : vi_suff,
+ 'pattern' : vi_patt, 'color' : vi_color
+ },
+ 'ndvi_aqua_sixteen_250':{'url':urlbase, 'folder':'MOLA/MYD13Q1.005',
+ 'res':250,'spec': vi_spec,'spec_qa': vi_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : vi_suff,
+ 'pattern' : vi_patt, 'color' : vi_color
+ },
+ 'ndvi_terra_sixteen_500':{'url':urlbase, 'folder':'MOLT/MOD13A2.005',
+ 'res':500,'spec': vi_spec,'spec_qa': vi_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : vi_suff,
+ 'pattern' : vi_patt, 'color' : vi_color
+ },
+ 'ndvi_aqua_sixteen_500':{'url':urlbase, 'folder':'MOLA/MYD13A2.005',
+ 'res':500,'spec': vi_spec,'spec_qa': vi_specqa,
+ 'spec_all' : lstvi_specall, 'suff' : vi_suff,
+ 'pattern' : vi_patt, 'color' : vi_color
+ }
}
+ snow = {'snow_terra_daily_500' : {'url' : usrsnow, 'folder' :
+ 'SAN/MOST/MOD10A1.005', 'res' : 500, 'spec' : snow_spec
+ ,'spec_qa': None},
+ 'snow_aqua_daily_500' : {'url' : usrsnow, 'folder' :
+ 'SAN/MOSA/MYD10A1.005', 'res' : 500, 'spec' : snow_spec
+ ,'spec_qa': None},
+ 'snow_terra_eight_500' : {'url' : usrsnow, 'folder' :
+ 'SAN/MOST/MOD10A2.005', 'res' : 500, 'spec' : snow_spec
+ ,'spec_qa': None},
+ 'snow_aqua_eight_500' : {'url' : usrsnow, 'folder' :
+ 'SAN/MOSA/MYD10A2.005', 'res' : 500, 'spec' : snow_spec
+ ,'spec_qa': None}
+ }
+ self.products = { }
+ self.products.update(lst)
+ self.products.update(vi)
+ self.products.update(snow)
+ self.products_swath = {
+ 'lst_terra_daily':{'url':urlbase,'folder':'MOLT/MOD11_L2.005',
+ 'spec': lstL2_spec}, 'lst_aqua_daily':{'url':urlbase,'folder':'MOLA/MYD11_L2.005',
+ 'spec': lstL2_spec}
+ }
def returned(self):
if self.products.keys().count(self.prod) == 1:
return self.products[self.prod]
@@ -83,6 +154,25 @@
grass.fatal(_("The code insert is not supported yet. Consider to ask " \
"on the grass-dev mailing list for future support"))
+ def color(self,code = None):
+ if code:
+ return self.fromcode(code)['color']
+ else:
+ return self.returned()['color']
+
+
+ def pattern(self,code = None):
+ if code:
+ return self.fromcode(code)['pattern']
+ else:
+ return self.returned()['pattern']
+
+ def suffix(self,code = None):
+ if code:
+ return self.fromcode(code)['suff']
+ else:
+ return self.returned()['suff']
+
def __str__(self):
prod = self.returned()
string = "url: " + prod['url'] + ", folder: " + prod['folder']
@@ -148,7 +238,10 @@
return self._outpar(0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, swath)
elif self.val == 'lcc':
SMajor = self._par('+a')
- SMinor = self._par('+b')
+ if self._par('+b'):
+ SMinor = self._par('+b')
+ elif self._par('+rf'):
+ SMinor = self._par('+rf')
STDPR1 = self._par('+lat_1')
STDPR2 = self._par('+lat_2')
CentMer = self._par('+lon_0')
@@ -159,7 +252,10 @@
CentLat, FE, FN, swath)
elif self.val == 'merc' or self.val == 'polar' or self.val == 'tmerc':
SMajor = self._par('+a')
- SMinor = self._par('+b')
+ if self._par('+b'):
+ SMinor = self._par('+b')
+ elif self._par('+rf'):
+ SMinor = self._par('+rf')
CentMer = self._par('+lon_0')
if self.val == 'tmerc':
Factor = self._par('+k_0')
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 2011-08-14 17:55:15 UTC (rev 47638)
+++ grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py 2011-08-15 07:56:40 UTC (rev 47639)
@@ -22,7 +22,7 @@
#############################################################################
#%module
-#% description: Download several tiles of MODIS products using pyModis
+#% description: Import single or more tiles of MODIS products using pyModis/MRT
#% keywords: raster
#% keywords: MODIS
#%end
@@ -36,13 +36,18 @@
#%end
#%flag
#% key: q
-#% description: Create and import also the QA layer
+#% description: Doesn't use the QA map
#%end
+#%flag
+#% key: r
+#% description: Doesn't resampling the output map
+#%end
#%option
#% key: mrtpath
#% type: string
#% key_desc: path
#% description: The full path to MRT directory
+#% gisprompt: old,dir,input
#% required: yes
#%end
#%option
@@ -50,6 +55,7 @@
#% type: string
#% key_desc: path
#% description: Path to single HDF file
+#% gisprompt: old,file,input
#% required: no
#%end
#%option
@@ -57,8 +63,13 @@
#% type: string
#% key_desc: file
#% description: Path to a file with a list of HDF files
+#% gisprompt: old,file,input
#% required: no
#%end
+#%option G_OPT_R_OUTPUT
+#% required : no
+#% guisection: Import
+#%end
#%option
#% key: resampl
#% type: string
@@ -75,16 +86,10 @@
#% description: A string to choose the subset of HDF file to use
#% required: no
#%end
-#%option
-#% key: folder
-#% type: string
-#% description: Folder where saves the data, full path
-#% answer: $HOME/.grass7/r.modis/import
-#% required: no
-#%end
import os, sys, string, glob, shutil
import grass.script as grass
+from datetime import date
# add the folder containing libraries to python path
libmodis = os.path.join(os.getenv('GISBASE'), 'etc', 'r.modis')
sys.path.append(libmodis)
@@ -95,7 +100,7 @@
except ImportError:
pass
-def list_files(opt, mosaik=False):
+def list_files(opt, mosaik = False):
"""Return a list of hdf files from the filelist on function single
Return a dictionary with a list o hdf file for each day on function mosaic
"""
@@ -126,7 +131,7 @@
basedir = os.path.split(filelist[0])[0]
return filelist, basedir
-def spectral(opts,code,q):
+def spectral(opts, code, q, m = False):
"""Return spectral string"""
# return the spectral set by the user
if opts['spectral'] != '':
@@ -134,16 +139,19 @@
# return the spectral by default
else:
prod = product().fromcode(code)
- if q:
- spectr = prod['spec_qa']
+ if m:
+ spectr = prod['spec_all']
+ elif q:
+ if prod['spec_qa']:
+ spectr = prod['spec_qa']
+ else:
+ spectr = prod['spec']
else:
spectr = prod['spec']
return spectr
-def confile(hdf,opts,q,mosaik=False):
+def confile(pm, opts, q, mosaik=False):
"""Create the configuration file for MRT software"""
- # create parseModis class and the parameter file
- pm = parseModis(hdf)
# return projection and datum
projObj = projection()
proj = projObj.returned()
@@ -152,24 +160,30 @@
zone = projObj.utmzone()
else:
zone = None
- cod = os.path.split(hdf)[1].split('.')[0]
+ cod = os.path.split(pm.hdfname)[1].split('.')[0]
if mosaik:
# if mosaic import all because the layers are choosen during mosaic
- spectr = "( 1 1 1 1 1 1 1 1 1 1 1 1 )"
+ spectr = spectral(opts, cod, q, True)
else:
- spectr = spectral(opts, cod,q)
+ spectr = spectral(opts, cod, q)
+ # out prefix
+ pref = prefix(opts)
# resampling
resampl = resampling(opts['resampl']).returned()
# projpar
projpar = projObj.return_params()
- return pm.confResample(spectr, None, None, dat, resampl, proj, zone, projpar)
+ return pm.confResample(spectr, None, pref, dat, resampl, proj, zone, projpar)
-def import_tif(hdf,basedir,rem,target=False):
+def prefix(options, name = False):
+ if options['output']:
+ return options['output']
+ else:
+ return None
+
+def import_tif(out, basedir, rem, target=False):
"""Import TIF files"""
- # prefix of HDF file
- prefix = os.path.split(hdf)[1].rstrip('.hdf')
# list of tif files
- tifiles = glob.glob1(basedir, prefix + "*.tif")
+ tifiles = glob.glob1(basedir, out + "*.tif")
if not tifiles:
grass.fatal(_('Error during the conversion'))
# check if is in latlong location to set flag l
@@ -177,6 +191,7 @@
f = "l"
else:
f = None
+ outfile = []
# for each file import it
for t in tifiles:
basename = os.path.splitext(t)[0]
@@ -184,15 +199,83 @@
try:
grass.run_command('r.in.gdal', input = name,
output = basename, flags = f, overwrite = True, quiet = True)
- if rem:
- os.remove(name)
- if target:
- shutil.move(name,target)
+ outfile.append(basename)
except:
grass.fatal(_('Error during import'))
- return 0
+ if rem:
+ os.remove(name)
+ if target != basedir:
+ shutil.move(name,target)
+ return outfile
-def single(options,remove,qa):
+def findfile(pref, suff):
+ if grass.find_file(pref + suff)['file']:
+ return grass.find_file(pref + suff)
+ else:
+ grass.fatal(_("Raster map <%s> not found") % pref + suff)
+
+def metadata(pars, mapp, coll):
+ meta = pars.metastring()
+ grass.run_command('r.support', map=mapp, hist=meta)
+ 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"))
+ if coll:
+ grass.run_command('r.colors', map=mapp, color = coll)
+
+def analize(pref, an, cod, parse):
+ prod = product().fromcode(cod)
+ if not prod['spec_qa']:
+ grass.warning(_("There is not QA file, analysis will be skipped"))
+ an = 'noqa'
+ pat = prod['pattern']
+ suf = prod['suff']
+ col = prod['color']
+ val = []
+ qa = []
+ for v,q in suf.iteritems():
+ val.append(findfile(pref,v))
+ qa.append(findfile(pref,q))
+
+ grass.run_command('g.region', rast = val[0]['fullname'])
+
+ if len(qa) != len(val):
+ grass.fatal(_("The number of QA and value maps is different, something wrong"))
+
+ for n in range(len(qa)):
+ valname = val[n]['name']
+ valfull = val[n]['fullname']
+ qaname = qa[n]['name']
+ qafull = qa[n]['fullname']
+ grass.run_command('r.null', map = valfull)
+ mapc = "%s.2 = (%s * 0.0200) - 273.15" % (valname, valfull)
+ grass.mapcalc(mapc)
+ if an == 'noqa':
+ grass.run_command('g.remove', rast = valfull)
+ grass.run_command('g.rename', rast= (valname + '.2', valfull))
+ #TODO check in modis.py to adjust the xml file of mosaic
+ #metadata(parse, valfull, col)
+ #metadata(parse, qafull, 'byr')
+
+ if an == 'all':
+ qamapcal = "%s.badpixels = if( " % qaname
+ for p in pat:
+ qamapcal += "%s%s || " % (qaname, p)
+ qamapcal.strip(' || ')
+ qamapcal += "1, null() )"
+ grass.mapcalc(qamapcal)
+ finalmap = "%s.3=if( isnull(%s.badpixels), %s.2, null())" % (
+ valname, qaname, valname)
+ grass.mapcalc(finalmap)
+ grass.run_command('g.remove', rast=(qaname + '.badpixels', valname,
+ valname + '.2'))
+ grass.run_command('g.rename', rast=(valname + '.3', valname))
+ #TODO check in modis.py to adjust the xml file of mosaic
+ #metadata(parse, valfull, col)
+ #metadata(parse, qafull, 'byr')
+
+def single(options,remove,an):
"""Convert the hdf file to tif and import it
"""
listfile, basedir = list_files(options)
@@ -201,27 +284,30 @@
# the full path to hdf file
hdf = os.path.join(basedir,i)
# create conf file fro mrt tools
- confname = confile(hdf,options,qa)
+ pm = parseModis(hdf)
+ confname = confile(pm,options,an)
# create convertModis class and convert it in tif file
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
- import_tif(i,basedir,remove)
+ maps_import = import_tif(output,basedir,remove)
+ if an:
+ cod = os.path.split(pm.hdfname)[1].split('.')[0]
+ analize(output, an, cod, pm)
os.remove(confname)
-
return grass.message(_('All files imported correctly'))
-#def createXMLmosaic(listfiles):
- #modis.
-
-def mosaic(options,remove,qa):
+def mosaic(options,remove,an):
"""Create a daily mosaic of hdf files convert to tif and import it
"""
dictfile, targetdir = list_files(options,True)
# for each day
for dat, listfiles in dictfile.iteritems():
# create the file with the list of name
- tempfile = open(grass.tempfile(),'w')
+ tempfile = open(os.path.join(targetdir,str(os.getpid())),'w')
tempfile.writelines(listfiles)
tempfile.close()
# basedir of tempfile, where hdf files are write
@@ -229,32 +315,38 @@
outname = "%s.%s.mosaic" % (listfiles[0].split('/')[-1].split('.')[0],
listfiles[0].split('/')[-1].split('.')[1])
# return the spectral subset in according mrtmosaic tool format
- spectr = spectral(options,listfiles[0].split('/')[-1].split('.')[0],qa)
+ spectr = spectral(options,listfiles[0].split('/')[-1].split('.')[0],an, True)
spectr = spectr.lstrip('( ').rstrip(' )')
# create mosaic
cm = createMosaic(tempfile.name,outname,options['mrtpath'],spectr)
cm.run()
# list of hdf files
- hdfdir = os.path.split(tempfile.name)[0]
- hdfiles = glob.glob1(hdfdir, outname + "*.hdf")
+ hdfiles = glob.glob1(basedir, outname + "*.hdf")
for i in hdfiles:
# the full path to hdf file
- hdf = os.path.join(hdfdir,i)
+ hdf = os.path.join(basedir,i)
# create conf file fro mrt tools
- confname = confile(hdf,options,qa,True)
+ pm = parseModis(hdf)
+ confname = confile(pm,options,an)
# create convertModis class and convert it in tif file
execmodis = convertModis(hdf,confname,options['mrtpath'])
execmodis.run()
# remove hdf
if remove:
# import tif files
- import_tif(i,basedir,remove)
+ import_tif(outname,basedir,remove)
+ if an:
+ cod = os.path.split(pm.hdfname)[1].split('.')[0]
+ analize(outname, an, cod, pm)
os.remove(hdf)
os.remove(hdf + '.xml')
# or move the hdf and hdf.xml to the dir where are the original files
else:
# import tif files
- import_tif(i,basedir,remove,targetdir)
+ import_tif(outname,basedir,remove,targetdir)
+ if an:
+ cod = os.path.split(pm.hdfname)[1].split('.')[0]
+ analize(outname, an, cod, pm)
try:
shutil.move(hdf,targetdir)
shutil.move(hdf + '.xml',targetdir)
@@ -288,24 +380,26 @@
if version['version'].find('7.') == -1:
grass.fatal(_('You are not in GRASS GIS version 7'))
return 0
- # check if also qa layer are converted
- if flags['q']:
- qa_layer = True
- else:
- qa_layer = False
# check if remove the files or not
if flags['t']:
remove = False
else:
remove = True
+ # check if do check quality, resampling and setting of color
+ if flags['r']:
+ analize = None
+ if flags['q']:
+ analize = 'noqa'
+ else:
+ analize = 'all'
# check if import simple file or mosaic
if flags['m'] and options['dns'] != '':
grass.fatal(_('It is not possible to create a mosaic with a single HDF file'))
return 0
elif flags['m']:
- mosaic(options,remove,qa_layer)
+ mosaic(options,remove,analize)
else:
- single(options,remove,qa_layer)
+ single(options,remove,analize)
if __name__ == "__main__":
options, flags = grass.parser()
More information about the grass-commit
mailing list