[GRASS-SVN] r46937 - in grass-addons/grass7/raster/r.modis: . libmodis r.modis.download r.modis.import

svn_grass at osgeo.org svn_grass at osgeo.org
Sun Jul 3 12:52:00 EDT 2011


Author: lucadelu
Date: 2011-07-03 09:52:00 -0700 (Sun, 03 Jul 2011)
New Revision: 46937

Added:
   grass-addons/grass7/raster/r.modis/r.modis.import/
   grass-addons/grass7/raster/r.modis/r.modis.import/Makefile
   grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.html
   grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py
Modified:
   grass-addons/grass7/raster/r.modis/Makefile
   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.download/r.modis.download.html
   grass-addons/grass7/raster/r.modis/r.modis.html
Log:
first version of r.modis.import

Modified: grass-addons/grass7/raster/r.modis/Makefile
===================================================================
--- grass-addons/grass7/raster/r.modis/Makefile	2011-07-03 16:35:23 UTC (rev 46936)
+++ grass-addons/grass7/raster/r.modis/Makefile	2011-07-03 16:52:00 UTC (rev 46937)
@@ -4,6 +4,7 @@
 
 SUBDIRS = \
           r.modis.download \
+          r.modis.import \
 	  libmodis
 
 include $(MODULE_TOPDIR)/include/Make/Dir.make

Modified: grass-addons/grass7/raster/r.modis/libmodis/modis.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/modis.py	2011-07-03 16:35:23 UTC (rev 46936)
+++ grass-addons/grass7/raster/r.modis/libmodis/modis.py	2011-07-03 16:52:00 UTC (rev 46937)
@@ -407,10 +407,18 @@
        filename = the name of MODIS hdf file
     """
     from xml.etree import ElementTree
-    # hdf name
-    self.hdfname = filename
-    # xml hdf name
-    self.xmlname = self.hdfname + '.xml'
+    if os.path.exists(filename):
+      # hdf name
+      self.hdfname = filename
+    else:
+      raise IOError('%s not exists' % self.hdfname)
+
+    if os.path.exists(self.hdfname + '.xml'):
+      # xml hdf name
+      self.xmlname = self.hdfname + '.xml'
+    else:
+      raise IOError('%s not exists' % self.hdfname + '.xml')
+
     # tif name for the output file for resample MRT software
     self.tifname = self.hdfname.replace('.hdf','.tif')
     with open(self.xmlname) as f:
@@ -454,22 +462,22 @@
   def retGranuleUR(self):
     """Return the GranuleUR element"""
     self.getGranule()
-    return self.granule.find('GranuleUR')
+    return self.granule.find('GranuleUR').text
 
   def retDbID(self):
     """Return the DbID element"""
     self.getGranule()
-    return self.granule.find('DbID')
+    return self.granule.find('DbID').text
 
   def retInsertTime(self):
     """Return the DbID element"""
     self.getGranule()
-    return self.granule.find('InsertTime')
+    return self.granule.find('InsertTime').text
 
   def retLastUpdate(self):
     """Return the DbID element"""
     self.getGranule()
-    return self.granule.find('LastUpdate')
+    return self.granule.find('LastUpdate').text
 
   def retCollectionMetaData(self):
     """Return the CollectionMetaData element"""
@@ -539,13 +547,14 @@
     value = {}
     self.getGranule()
     mes = self.granule.find('MeasuredParameter')
-    value['ParameterName'] = mes.find('ParameterName').text
-    meStat = mes.find('MeasuredParameterContainer').find('QAStats')
+    mespc = mes.find('MeasuredParameterContainer')
+    value['ParameterName'] = mespc.find('ParameterName').text
+    meStat = mespc.find('QAStats')
     qastat = {}
     for i in meStat.getiterator():
       qastat[i.tag] = i.text
     value['QAStats'] = qastat
-    meFlag = mes.find('MeasuredParameterContainer').find('QAFlags')
+    meFlag = mespc.find('QAFlags')
     flagstat = {}
     for i in meStat.getiterator():
       flagstat[i.tag] = i.text
@@ -556,8 +565,9 @@
     """Return the platform values inside a dictionary."""
     value = {}
     self.getGranule()
-    value['PlatformShortName'] = self.granule.find('PlatformShortName').text
-    instr = self.granule.find('Instrument')
+    plat = self.granule.find('Platform')
+    value['PlatformShortName'] = plat.find('PlatformShortName').text
+    instr = plat.find('Instrument')
     value['InstrumentShortName'] = instr.find('InstrumentShortName').text
     sensor = instr.find('Sensor')
     value['SensorShortName'] = sensor.find('SensorShortName').text
@@ -734,8 +744,184 @@
     else:
       subprocess.call([execut,'-p',self.conf])
     return "The hdf file %s was converted" % self.name
-      
-#class createMosaic:
-  #def __init__(self
-              #listfile,
-              #mrtpath)
+
+class createMosaic:
+  """A class to convert a mosaic of different modis tiles"""
+  def __init__(self,
+              listfile,
+              outprefix,
+              mrtpath,
+              subset = False):
+    # check if the hdf file exists
+    if os.path.exists(listfile):
+      self.basepath = os.path.split(listfile)[0]
+      self.listfiles = listfile
+      self.HDFfiles = open(listfile).readlines()
+    else:
+      raise IOError('%s not exists' % hdfname)
+    # check if mrtpath and subdirectories exists and set environment variables
+    if os.path.exists(mrtpath):
+      if os.path.exists(os.path.join(mrtpath,'bin')):
+        self.mrtpathbin = os.path.join(mrtpath,'bin')
+        os.environ['PATH'] = "%s:%s" % (os.environ['PATH'],os.path.join(mrtpath,
+                                                                        'data'))
+      else:
+        raise IOError('The path %s not exists' % os.path.join(mrtpath,'bin'))
+      if os.path.exists(os.path.join(mrtpath,'data')):
+        self.mrtpathdata = os.path.join(mrtpath,'data')
+        os.environ['MRTDATADIR'] = os.path.join(mrtpath,'data')
+      else:
+        raise IOError('The path %s not exists' % os.path.join(mrtpath,'data'))
+    else:
+      raise IOError('The path %s not exists' % mrtpath)
+    self.out = os.path.join(self.basepath, outprefix + '.hdf')
+    self.outxml = os.path.join(self.basepath, self.out + '.xml')
+    self.subset = subset
+
+  def boundaries(self):
+    """Return the max extend for the mosaic"""
+    pm = parseModis(os.path.join(self.basepath,self.HDFfiles[0].strip()))
+    boundary = pm.retBoundary()
+    for i in range(1,len(self.HDFfiles)):
+      pm = parseModis(os.path.join(self.basepath,self.HDFfiles[i].strip()))
+      bound = pm.retBoundary()
+      if bound['min_lat'] < boundary['min_lat']:
+        boundary['min_lat'] = bound['min_lat']
+      if bound['min_lon'] < boundary['min_lon']:
+        boundary['min_lon'] = bound['min_lon']
+      if bound['max_lat'] > boundary['max_lat']:
+        boundary['max_lat'] = bound['max_lat']
+      if bound['max_lon'] > boundary['max_lon']:
+        boundary['max_lon'] = bound['max_lon']
+    return boundary
+
+  def write_mosaic_xml(self):
+    """Return a xml file for a mosaic"""
+    from xml.etree import ElementTree
+
+    def cicle_values(ele,values):
+      """Function to add values from a dictionary"""
+      for k,v in values.iteritems():
+        elem = ElementTree.SubElement(ele,k)
+        elem.text = v
+
+    # take the data similar for different files from the first
+    pm = parseModis(os.path.join(self.basepath,self.HDFfiles[0].strip()))
+    # the root element
+    granule = ElementTree.Element('GranuleMetaDataFile')
+    # add DTDVersion
+    dtd = ElementTree.SubElement(granule,'DTDVersion')
+    dtd.text = pm.retDTD()    
+    # add DataCenterId
+    dci = ElementTree.SubElement(granule,'DataCenterId')
+    dci.text = pm.retDataCenter()
+    # add GranuleURMetaData
+    gurmd = ElementTree.SubElement(granule,'GranuleURMetaData')
+
+    #gur = ElementTree.SubElement(gurmd,'GranuleUR')
+    #gur.text = pm.retGranuleUR()
+    # TODO ADD GranuleUR DbID InsertTime LastUpdate
+
+    # add CollectionMetaData
+    coll = pm.retCollectionMetaData()
+    cmd = ElementTree.SubElement(gurmd,'CollectionMetaData')
+    cicle_values(cmd,coll)
+
+    ## DA RIVEDERE DataFileContainer
+    #dataf = pm.retDataFiles()
+    #df = ElementTree.SubElement(gurmd,'CollectionMetaData')
+    #dfc = ElementTree.SubElement(df,'DataFileContainer')
+    #cicle_values(dfc,dataf)
+
+    # add PGEVersionClass
+    pgevc = ElementTree.SubElement(gurmd,'PGEVersionClass')
+    pgevc.text = pm.retPGEVersion()
+    # add RangeDateTime
+    datime = pm.retRangeTime()
+    rdt = ElementTree.SubElement(gurmd,'RangeDateTime')
+    cicle_values(rdt,datime)
+    # SpatialDomainContainer
+    maxbound = self.boundaries()
+    sdc = ElementTree.SubElement(gurmd,'SpatialDomainContainer')
+    hsdc = ElementTree.SubElement(sdc,'HorizontalSpatialDomainContainer')
+    gp = ElementTree.SubElement(hsdc,'GPolygon')
+    bound = ElementTree.SubElement(gp,'Boundary')
+    pt1 = ElementTree.SubElement(bound, 'Point')
+    pt1lon = ElementTree.SubElement(pt1, 'PointLongitude')
+    pt1lon.text = str(maxbound['min_lon'])
+    pt1lat = ElementTree.SubElement(pt1, 'PointLatitude')
+    pt1lat.text = str(maxbound['max_lat'])
+    pt2 = ElementTree.SubElement(bound, 'Point')
+    pt2lon = ElementTree.SubElement(pt2, 'PointLongitude')
+    pt2lon.text = str(maxbound['max_lon'])
+    pt2lat = ElementTree.SubElement(pt2, 'PointLatitude')
+    pt2lat.text = str(maxbound['max_lat'])
+    pt3 = ElementTree.SubElement(bound, 'Point')
+    pt3lon = ElementTree.SubElement(pt3, 'PointLongitude')
+    pt3lon.text = str(maxbound['min_lon'])
+    pt3lat = ElementTree.SubElement(pt3, 'PointLatitude')
+    pt3lat.text = str(maxbound['min_lat'])
+    pt4 = ElementTree.SubElement(bound, 'Point')
+    pt4lon = ElementTree.SubElement(pt4, 'PointLongitude')
+    pt4lon.text = str(maxbound['max_lon'])
+    pt4lat = ElementTree.SubElement(pt4, 'PointLatitude')
+    pt4lat.text = str(maxbound['min_lat'])
+    # add MeasuredParameter
+    mp = ElementTree.SubElement(gurmd,'MeasuredParameter')
+    mpc = ElementTree.SubElement(mp,'MeasuredParameterContainer')
+    pn = ElementTree.SubElement(mpc,'ParameterName')
+    measure = pm.retMeasure()
+    pn.text = measure['ParameterName']
+
+    # TODO ADD qstats qflags
+    # add Platform
+    platvalues = pm.retPlatform()
+    plat = ElementTree.SubElement(gurmd,'Platform')
+    psn = ElementTree.SubElement(plat,'PlatformShortName')
+    psn.text = platvalues['PlatformShortName']
+    ins = ElementTree.SubElement(plat,'Instrument')
+    isn = ElementTree.SubElement(ins,'InstrumentShortName')
+    isn.text = platvalues['InstrumentShortName']
+    sens = ElementTree.SubElement(ins,'Sensor')
+    ssn = ElementTree.SubElement(sens,'SensorShortName')
+    ssn.text = platvalues['SensorShortName']
+    # add PSAs
+    psas = ElementTree.SubElement(gurmd,'PSAs')
+    # TODO ADD all PSA
+    ig = ElementTree.SubElement(gurmd,'InputGranule')
+    for i in self.HDFfiles:
+      ip = ElementTree.SubElement(ig,'InputPointer')
+      ip.text = i
+    # TODO ADD BrowseProduct
+    output = open(self.outxml, 'w')
+    output.write('<?xml version="1.0" encoding="UTF-8"?>')
+    output.write('<!DOCTYPE GranuleMetaDataFile SYSTEM "http://ecsinfo.gsfc.nasa.gov/ECSInfo/ecsmetadata/dtds/DPL/ECS/ScienceGranuleMetadata.dtd">')
+    output.write(ElementTree.tostring(granule))
+    output.close()
+
+  def executable(self):
+    """Return the executable of mrtmosaic MRT software
+    """
+    if sys.platform.count('linux') != -1:
+      if os.path.exists(os.path.join(self.mrtpathbin,'mrtmosaic')):
+        return os.path.join(self.mrtpathbin,'mrtmosaic')
+    elif sys.platform.count('win32') != -1:
+      if os.path.exists(os.path.join(self.mrtpathbin,'mrtmosaic.exe')):
+        return os.path.join(self.mrtpath,'mrtmosaic.exe')
+
+  def run(self):
+    """Exect the mosaic process"""
+    import subprocess
+    execut = self.executable()
+    if not os.path.exists(execut):
+      raise IOError('The path %s not exists, could be an erroneus path or '\
+                    + 'software') % execut
+    else:
+      self.write_mosaic_xml()
+      if self.subset:
+        subprocess.call([execut,'-i',self.listfiles,'-o',self.out,'-s',self.subset], 
+                        stderr = subprocess.STDOUT)
+      else:
+        subprocess.call([execut,'-i',self.listfiles,'-o',self.out], stderr = 
+                        subprocess.STDOUT)
+    return "The mosaic file %s is created" % self.out

Modified: grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py	2011-07-03 16:35:23 UTC (rev 46936)
+++ grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py	2011-07-03 16:52:00 UTC (rev 46937)
@@ -51,13 +51,13 @@
 	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':'( 0 1 0 0 0 1 0 0 0 0 0 0 )'},
+          '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':'( 0 1 0 0 0 1 0 0 0 0 0 0 )'}, 
+          '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':'( 0 0 1 1 0 0 0 0 0 0 0 )'}
+          '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 )'}
         }
 
     def returned(self):
@@ -94,7 +94,7 @@
     another one. Not all projection systems are supported"""
     def __init__(self,value):
         self.proj = value
-        self.projections = {'latlong':'GEO', 'lcc':'LAMBERT CONFORMAL CONIC',
+        self.projections = {'ll':'GEO', 'lcc':'LAMBERT CONFORMAL CONIC',
              'merc':'MERCATOR', 'polar':'POLAR STEREOGRAPHIC', 'utm':'UTM', 
              'tmerc':'TRANSVERSE MERCATOR'}
 

Modified: grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.html
===================================================================
--- grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.html	2011-07-03 16:35:23 UTC (rev 46936)
+++ grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.html	2011-07-03 16:52:00 UTC (rev 46937)
@@ -45,7 +45,7 @@
 <h2>SEE ALSO</h2>
 <ul>
   <li> <a href="r.modis.html">r.modis</a></li>
-<!--   <li> <a href="r.modis.import.html">r.modis.import</a></li> -->
+  <li> <a href="r.modis.import.html">r.modis.import</a></li>
 <!--  <li> <a href="r.modis.process.html">r.modis.process</a>: Calculates range of patch area size on a raster map</li>-->
 </ul>
  

Modified: grass-addons/grass7/raster/r.modis/r.modis.html
===================================================================
--- grass-addons/grass7/raster/r.modis/r.modis.html	2011-07-03 16:35:23 UTC (rev 46936)
+++ grass-addons/grass7/raster/r.modis/r.modis.html	2011-07-03 16:52:00 UTC (rev 46937)
@@ -38,8 +38,8 @@
 <h2>SEE ALSO</h2>
 <ul>
   <li> <a href="r.modis.download.html">r.modis.download</a>: Download MODIS HDF products from NASA FTP</li>
-<!--   <li> <a href="r.modis.import.html">r.modis.import</a>: Import Level 3 MODIS products as 
-single image or daily mosaik into GRASS GIS</li> -->
+  <li> <a href="r.modis.import.html">r.modis.import</a>: Import Level 3 MODIS products as 
+single image or daily mosaik into GRASS GIS</li>
 <!--  <li> <a href="r.modis.process.html">r.modis.process</a>: Calculates range
 of patch area size on a raster map</li>-->
 </ul>

Added: grass-addons/grass7/raster/r.modis/r.modis.import/Makefile
===================================================================
--- grass-addons/grass7/raster/r.modis/r.modis.import/Makefile	                        (rev 0)
+++ grass-addons/grass7/raster/r.modis/r.modis.import/Makefile	2011-07-03 16:52:00 UTC (rev 46937)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../../..
+
+PGM = r.modis.import
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
+
+# $(ETC)/r.modis/r.modis.download/%: % | $(ETC)/r.modis/r.modis.download
+# 	$(INSTALL) $< $@
+# 
+# $(ETC)/r.modis/r.modis.download:
+# 	$(MKDIR) $@
\ No newline at end of file

Added: grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.html
===================================================================
--- grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.html	                        (rev 0)
+++ grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.html	2011-07-03 16:52:00 UTC (rev 46937)
@@ -0,0 +1,48 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.modis.import</em> import into GRASS GIS Level 3 MODIS product.
+
+<h2>NOTE</h2>
+
+The input file is the list of full path to HDF files, one for line.
+The input file have to be inside the folder where the HDF files are.<br>
+The mrtpath option is the path to the main folder of Modis Reprojection Tools; it contain bin and data folder, this two folders are essential for a successful result.<br>
+It is not possible to use flag "q" and option "spectral" together.
+<br>
+<br>
+
+<b>Warning</b>: 
+<ul>
+  <li>Currently it's tested and It work correctly only with Lat/Long</li>
+  <li>Currently the mosaic flag it is not working</li>
+</ul>
+
+<h2>EXAMPLES</h2>
+
+Import file with the default subset without QA layer
+<div class="code"><pre>
+    r.modis.import filename=/path/to/listfile mrtpath=/path/to/mrt
+</pre></div>
+
+Import file with the default subset with QA layer
+<div class="code"><pre>
+    r.modis.import -q filename=/path/to/listfile mrtpath=/path/to/mrt
+</pre></div>
+
+Import file with the your own subset
+<div class="code"><pre>
+    r.modis.import filename=/path/to/listfile mrtpath=/path/to/mrt spectral="( 1 )"
+</pre></div>
+
+<h2>SEE ALSO</h2>
+<ul>
+  <li> <a href="r.modis.html">r.modis</a></li>
+  <li> <a href="r.modis.download.html">r.modis.download</a></li>
+<!--  <li> <a href="r.modis.process.html">r.modis.process</a>: Calculates range of patch area size on a raster map</li>-->
+</ul>
+ 
+<h2>AUTHOR</h2>
+
+Luca Delucchi, Google Summer of Code 2011
+
+<p><i>Last changed: $Date: 2011-06-18 14:44:59 +0200 (Mon, 18 June 2011) $</i>
\ No newline at end of file

Added: 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	                        (rev 0)
+++ grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py	2011-07-03 16:52:00 UTC (rev 46937)
@@ -0,0 +1,302 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+############################################################################
+#
+# MODULE:        r.modis.import
+# AUTHOR(S):     Luca Delucchi
+# PURPOSE:       r.modis.import is an internafe to pyModis to download 
+#                several tiles of MODIS produts from NASA ftp
+#
+# COPYRIGHT:        (C) 2011 by Luca Delucchi
+#
+#                This program is free software under the GNU General Public
+#                License (>=v2). Read the file COPYING that comes with GRASS
+#                for details.
+#
+#############################################################################
+#
+# REQUIREMENTS:
+#      -  MRT tools
+#
+#############################################################################
+
+#%module
+#% description: Download several tiles of MODIS products using pyModis
+#% keywords: raster
+#% keywords: MODIS
+#%end
+#%flag
+#% key: m
+#% description: Create a mosaic for each day
+#%end
+#%flag
+#% key: t
+#% description: Preserve temporary files (TIF and HDF mosaic)
+#%end
+#%flag
+#% key: q
+#% description: Create and import also the QA layer
+#%end
+#%option
+#% key: mrtpath
+#% type: string
+#% key_desc: path
+#% description: The full path to MRT directory
+#% required: yes
+#%end
+#%option
+#% key: dns
+#% type: string
+#% key_desc: path
+#% description: Path to single HDF file
+#% required: no
+#%end
+#%option
+#% key: files
+#% type: string
+#% key_desc: file
+#% description: Path to a file with a list of HDF files 
+#% required: no
+#%end
+#%option
+#% key: resampl
+#% type: string
+#% key_desc: resampling
+#% description: Code of resampling type
+#% options: NN, BI, CC, NONE
+#% answer: NN
+#% required: no
+#%end
+#%option
+#% key: spectral
+#% type: string
+#% key_desc: spectral subset
+#% 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
+# add the folder containing libraries to python path
+libmodis = os.path.join(os.getenv('GISBASE'), 'etc', 'r.modis')
+sys.path.append(libmodis)
+# try to import pymodis (modis) and some class for r.modis.download
+try:
+    from rmodislib import resampling, product, get_proj, projection, datum
+    from modis  import parseModis, convertModis, createMosaic
+except ImportError:
+    pass
+
+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
+    """
+    # read the file with the list of HDF
+    if opt['files'] != '':
+        listoffile = open(opt['files'],'r')
+        basedir = os.path.split(listoffile.name)[0]
+        # if mosaic create a dictionary
+        if mosaik:
+            filelist = {}
+        # if not mosaic create a list
+        else:
+            filelist = []
+        # append hdf files
+        for line in listoffile:
+            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:
+                day = line.split('/')[-1].split('.')[1]
+                if filelist.has_key(day):
+                    filelist[day].append(line)
+                else:
+                    filelist[day] = [line]
+    # create a list for each file
+    elif options['dns'] != '':
+        filelist = [options['dns']]
+        basedir = os.path.split(filelist[0])[0]
+    return filelist, basedir
+
+def spectral(opts,code,q):
+    """Return spectral string"""
+    # return the spectral set by the user
+    if opts['spectral'] != '':
+        spectr = opts['spectral']
+    # return the spectral by default
+    else:
+        prod = product().fromcode(code)
+        if q:
+            spectr = prod['spec_qa']
+        else:
+            spectr = prod['spec']
+    return spectr
+
+def confile(hdf,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 
+    inproj = get_proj()
+    proj = projection(inproj['proj']).returned()
+    if proj != 'GEO':
+        grass.fatal(_('Projection different to Lat/Long is not working yet. Sorry'))
+        return 0
+    dat = datum(inproj['datum']).returned()
+    if proj == 'UTM':
+        zone = prod['zone']
+    else:
+        zone = None
+    cod = os.path.split(hdf)[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 )"
+    else:
+        spectr = spectral(opts, cod,q)
+    # resampling
+    resampl = resampling(opts['resampl']).returned()
+    # projpar TO DO
+    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 )'
+    return pm.confResample(spectr, None, None, dat, resampl, proj, zone, projpar)
+
+def import_tif(hdf,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")
+    # for each file import it
+    for t in tifiles:
+        basename = os.path.splitext(t)[0]
+        name = os.path.join(basedir,t)
+        grass.run_command('r.in.gdal', input = name, 
+        output = basename, flags = "l", overwrite = True, quiet = True)
+        if rem:
+            os.remove(name)
+        if target:
+             shutil.move(name,target)
+    return 0
+
+def single(options,remove,qa):
+    """Convert the hdf file to tif and import it
+    """
+    listfile, basedir = list_files(options)
+    # for each file
+    for i in listfile:
+        # the full path to hdf file
+        hdf = os.path.join(basedir,i)
+        # create conf file fro mrt tools
+        confname = confile(hdf,options,qa)
+        # create convertModis class and convert it in tif file
+        execmodis = convertModis(hdf,confname,options['mrtpath'])
+        execmodis.run()
+        # import tif files
+        import_tif(i,basedir,remove)
+        os.remove(confname)
+
+    return grass.message(_('All files imported correctly'))
+
+def mosaic(options,remove,qa):
+    """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.writelines(listfiles)
+        tempfile.close()
+        # basedir of tempfile, where hdf files are write
+        basedir = os.path.split(tempfile.name)[0]
+        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 = 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")
+        for i in hdfiles:
+            # the full path to hdf file
+            hdf = os.path.join(hdfdir,i)
+            # create conf file fro mrt tools
+            confname = confile(hdf,options,qa,True)
+            # 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)
+                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)
+                try: 
+                    shutil.move(hdf,targetdir)
+                    shutil.move(hdf + '.xml',targetdir)
+                except:
+                    pass
+            # remove the conf file
+            os.remove(confname)
+        grass.try_remove(tempfile.name)
+    return grass.message(_('All files imported correctly'))
+
+def main():
+    # check if you are in GRASS
+    gisbase = os.getenv('GISBASE')
+    if not gisbase:
+        grass.fatal(_('$GISBASE not defined'))
+        return 0
+    # return an error if q and spectral are set
+    if flags['q'] and options['spectral'] != '':
+        grass.fatal(_('It is not possible set flag "q" and "spectral" option'))
+        return 0
+    # return an error if both dns and files option are set or not
+    if options['dns'] == '' and options['files'] == '':
+        grass.fatal(_('You have to choese one of "dns" or "spectral" option'))
+        return 0
+    elif options['dns'] != '' and options['files'] != '':
+        grass.fatal(_('It is not possible set  "dns" and "spectral" option together'))
+        return 0
+    # check the version
+    version = grass.core.version()
+    # this is would be set automatically
+    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 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)
+    else:
+        single(options,remove,qa_layer)
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    sys.exit(main())



More information about the grass-commit mailing list