[GRASS-SVN] r66526 - in grass-addons/grass7/raster/r.modis: libmodis r.modis.download r.modis.import
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Oct 19 06:43:18 PDT 2015
Author: lucadelu
Date: 2015-10-19 06:43:18 -0700 (Mon, 19 Oct 2015)
New Revision: 66526
Added:
grass-addons/grass7/raster/r.modis/libmodis/convertmodis_gdal.py
Modified:
grass-addons/grass7/raster/r.modis/libmodis/Makefile
grass-addons/grass7/raster/r.modis/libmodis/downmodis.py
grass-addons/grass7/raster/r.modis/libmodis/parsemodis.py
grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.py
grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py
Log:
r.modis: several updates, r.modis.import to be tested
Modified: grass-addons/grass7/raster/r.modis/libmodis/Makefile
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/Makefile 2015-10-19 01:03:34 UTC (rev 66525)
+++ grass-addons/grass7/raster/r.modis/libmodis/Makefile 2015-10-19 13:43:18 UTC (rev 66526)
@@ -3,7 +3,7 @@
include $(MODULE_TOPDIR)/include/Make/Other.make
include $(MODULE_TOPDIR)/include/Make/Python.make
-MODULES = downmodis rmodislib convertmodis parsemodis
+MODULES = downmodis rmodislib convertmodis parsemodis convertmodis_gdal
ETCDIR = $(ETC)/r.modis
Added: grass-addons/grass7/raster/r.modis/libmodis/convertmodis_gdal.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/convertmodis_gdal.py (rev 0)
+++ grass-addons/grass7/raster/r.modis/libmodis/convertmodis_gdal.py 2015-10-19 13:43:18 UTC (rev 66526)
@@ -0,0 +1,710 @@
+#!/usr/bin/env python
+# class to convert/process modis data using gdal
+#
+# Some code is coming from GDAL/OGR Test Suite
+# Copyright (c) 2008, Andrey Kiselev <dron16 at ak4719.spb.edu>
+# Copyright (c) 2008-2014, Even Rouault <even dot rouault at mines-paris dot org>
+# http://trac.osgeo.org/gdal/browser/trunk/autotest/alg/warp.py#L782
+#
+# (c) Copyright Luca Delucchi 2014
+# Authors: Luca Delucchi
+# Email: luca dot delucchi at iasma dot it
+#
+##################################################################
+#
+# This MODIS Python class is licensed under the terms of GNU GPL 2.
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as
+# published by the Free Software Foundation; either version 2 of
+# the License, or (at your option) any later version.
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU General Public License for more details.
+#
+##################################################################
+"""Convert MODIS HDF file using GDAL Python bindings. It can create GeoTiff
+file (or other GDAL supported formats) or HDF mosaic file for several tiles.
+
+Classes:
+
+* :class:`file_info`
+* :class:`createMosaicGDAL`
+* :class:`convertModisGDAL`
+
+Functions:
+
+* :func:`getResampling`
+* :func:`raster_copy`
+* :func:`raster_copy_with_nodata`
+
+"""
+from __future__ import print_function
+from types import ListType, StringType
+try:
+ import osgeo.gdal as gdal
+except ImportError:
+ try:
+ import gdal
+ except ImportError:
+ raise ImportError('Python GDAL library not found, please install '
+ 'python-gdal')
+
+try:
+ import osgeo.osr as osr
+except ImportError:
+ try:
+ import osr
+ except ImportError:
+ raise ImportError('Python GDAL library not found, please install '
+ 'python-gdal')
+
+
+RESAM_GDAL = ['AVERAGE', 'BILINEAR', 'CUBIC', 'CUBIC_SPLINE', 'LANCZOS',
+ 'MODE', 'NEAREST_NEIGHBOR']
+SINU_WKT = 'PROJCS["Sinusoidal_Sanson_Flamsteed",GEOGCS["GCS_Unknown",' \
+ 'DATUM["D_unknown",SPHEROID["Unknown",6371007.181,"inf"]],' \
+ 'PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]' \
+ ',PROJECTION["Sinusoidal"],PARAMETER["central_meridian",0],' \
+ 'PARAMETER["false_easting",0],PARAMETER["false_northing",0]' \
+ ',UNIT["Meter",1]]'
+
+
+def getResampling(res):
+ """Return the GDAL resampling method
+
+ :param str res: the string of resampling method
+ """
+ if res == 'AVERAGE':
+ return gdal.GRA_Average
+ elif res == 'BILINEAR' or res == 'BICUBIC':
+ return gdal.GRA_Bilinear
+ elif res == 'LANCZOS':
+ return gdal.GRA_Lanczos
+ elif res == 'MODE':
+ return gdal.GRA_Mode
+ elif res == 'NEAREST_NEIGHBOR':
+ return gdal.GRA_NearestNeighbour
+ elif res == 'CUBIC_CONVOLUTION' or res == 'CUBIC':
+ return gdal.GRA_Cubic
+ elif res == 'CUBIC_SPLINE':
+ return gdal.GRA_CubicSpline
+
+
+class convertModisGDAL:
+ """A class to convert modis data from hdf to GDAL formats using GDAL
+
+ :param str hdfname: name of input data
+ :param str prefix: prefix for output data
+ :param str subset: the subset to consider
+ :param int res: output resolution
+ :param str outformat: output format, it is possible to use all the
+ supported GDAL format
+ :param int epsg: the EPSG code for the reprojection of output file
+ :param str wkt: the WKT string for the reprojection of output file
+ :param str resampl: the resampling method to use
+ :param bool vrt: True to read GDAL VRT file created with
+ createMosaicGDAL
+ """
+ def __init__(self, hdfname, prefix, subset, res, outformat="GTiff",
+ epsg=None, wkt=None, resampl='NEAREST_NEIGHBOR', vrt=False):
+ """Function for the initialize the object"""
+ # Open source dataset
+ self.in_name = hdfname
+ self.src_ds = gdal.Open(self.in_name)
+ self.layers = self.src_ds.GetSubDatasets()
+ self.output_pref = prefix
+ self.resolution = res
+ if epsg:
+ self.dst_srs = osr.SpatialReference()
+ self.dst_srs.ImportFromEPSG(int(epsg))
+ self.dst_wkt = self.dst_srs.ExportToWkt()
+ elif wkt:
+ try:
+ f = open(wkt)
+ self.dst_wkt = f.read()
+ f.close()
+ except:
+ self.dst_wkt = wkt
+ else:
+ raise Exception('You have to set one of the following option: '
+ '"epsg", "wkt"')
+ # error threshold the same value as gdalwarp
+ self.error_threshold = 0.125
+ self.resampling = getResampling(resampl)
+ if type(subset) == ListType:
+ self.subset = subset
+ elif type(subset) == StringType:
+ self.subset = subset.replace('(', '').replace(')',
+ '').strip().split()
+ else:
+ raise Exception('Type for subset parameter not supported')
+ self.driver = gdal.GetDriverByName(outformat)
+ self.vrt = vrt
+ if self.driver is None:
+ raise Exception('Format driver %s not found, pick a supported '
+ 'driver.' % outformat)
+
+ def _boundingBox(self, src):
+ """Obtain the bounding box of raster in the new coordinate system
+
+ :param src: a GDAL dataset object
+
+ :return: a bounding box value in lists
+ """
+ src_gtrn = src.GetGeoTransform(can_return_null=True)
+
+ src_bbox_cells = ((0., 0.), (0, src.RasterYSize), (src.RasterXSize, 0),
+ (src.RasterXSize, src.RasterYSize))
+
+ geo_pts_x = []
+ geo_pts_y = []
+ for x, y in src_bbox_cells:
+ x2 = src_gtrn[0] + src_gtrn[1] * x + src_gtrn[2] * y
+ y2 = src_gtrn[3] + src_gtrn[4] * x + src_gtrn[5] * y
+ geo_pts_x.append(x2)
+ geo_pts_y.append(y2)
+ return ((min(geo_pts_x), min(geo_pts_y)), (max(geo_pts_x),
+ max(geo_pts_y)))
+
+ def _calculateRes(self, minn, maxx, res):
+ """Calculate the number of pixel from extent and resolution
+
+ :param float minn: minimum value of extent
+ :param float maxx: maximum value of extent
+ :param int res: resolution of output raster
+
+ :return: integer number with the number of pixels
+ """
+ return int(round((maxx - minn) / res))
+
+ def _createWarped(self, raster):
+ """Create a warped VRT file to fetch default values for target raster
+ dimensions and geotransform
+
+ :param str raster: the name of raster, for HDF have to be one subset
+ """
+ src = gdal.Open(raster)
+ tmp_ds = gdal.AutoCreateWarpedVRT(src, SINU_WKT,
+ self.dst_wkt, self.resampling,
+ self.error_threshold)
+
+ if not self.resolution:
+ self.dst_xsize = tmp_ds.RasterXSize
+ self.dst_ysize = tmp_ds.RasterYSize
+ self.dst_gt = tmp_ds.GetGeoTransform()
+ else:
+ bbox = self._boundingBox(tmp_ds)
+ self.dst_xsize = self._calculateRes(bbox[0][0], bbox[1][0],
+ self.resolution)
+ self.dst_ysize = self._calculateRes(bbox[0][1], bbox[1][1],
+ self.resolution)
+ if self.dst_xsize == 0:
+ raise Exception('Invalid number of pixel 0 for X size. The '
+ 'problem could be in an invalid value of '
+ 'resolution')
+ return 0
+ elif self.dst_ysize == 0:
+ raise Exception('Invalid number of pixel 0 for Y size. The '
+ 'problem could be in an invalid value of '
+ 'resolution')
+ return 0
+ self.dst_gt = [bbox[0][0], self.resolution, 0.0, bbox[1][1], 0.0,
+ -self.resolution]
+ tmp_ds = None
+ src = None
+ return 0
+
+ def _progressCallback(self, pct, message, user_data):
+ """For the progress status"""
+ return 1 # 1 to continue, 0 to stop
+
+ def _reprojectOne(self, l):
+ """Reproject a single subset of MODIS product
+
+ l = complete name of input dataset
+ """
+ l_src_ds = gdal.Open(l)
+ meta = l_src_ds.GetMetadata()
+ band = l_src_ds.GetRasterBand(1)
+ if '_FillValue' in meta.keys():
+ fill_value = meta['_FillValue']
+ elif band.GetNoDataValue():
+ fill_value = band.GetNoDataValue()
+ else:
+ fill_value = None
+ datatype = band.DataType
+ try:
+ l_name = l.split(':')[-1]
+ out_name = "{pref}_{lay}.tif".format(pref=self.output_pref,
+ lay=l_name)
+ except:
+ out_name = "{pref}.tif".format(pref=self.output_pref)
+ if self.vrt:
+ out_name = "{pref}.tif".format(pref=self.output_pref)
+ try:
+ dst_ds = self.driver.Create(out_name, self.dst_xsize,
+ self.dst_ysize, 1, datatype)
+ except:
+ raise Exception('Not possibile to create dataset %s' % out_name)
+ return 0
+ dst_ds.SetProjection(self.dst_wkt)
+ dst_ds.SetGeoTransform(self.dst_gt)
+ if fill_value:
+ dst_ds.GetRasterBand(1).SetNoDataValue(float(fill_value))
+ dst_ds.GetRasterBand(1).Fill(float(fill_value))
+ cbk = self._progressCallback
+ # value for last parameter of above self._progressCallback
+ cbk_user_data = None
+ try:
+ gdal.ReprojectImage(l_src_ds, dst_ds, SINU_WKT, self.dst_wkt,
+ self.resampling, 0, self.error_threshold, cbk,
+ cbk_user_data)
+ print("Layer {name} reprojected".format(name=l))
+ except:
+ raise Exception('Not possibile to reproject dataset '
+ '{name}'.format(name=l))
+ return 0
+ dst_ds.SetMetadata(meta)
+ dst_ds = None
+ l_src_ds = None
+ return 0
+
+ def run_vrt_separated(self):
+ """Reproject VRT created by createMosaicGDAL, function write_vrt with
+ sepatated=True
+ """
+ self._createWarped(self.in_name)
+ self._reprojectOne(self.in_name)
+ print("Dataset '{name}' reprojected".format(name=self.in_name))
+
+ def run(self):
+ """Reproject all the subset of choosen layer"""
+ if self.vrt:
+ self.run_vrt_separated()
+ return
+ else:
+ self._createWarped(self.layers[0][0])
+ n = 0
+ for i in self.subset:
+ if str(i) == '1':
+ self._reprojectOne(self.layers[n][0])
+ n = n + 1
+ print("All layer for dataset '{name}' "
+ "reprojected".format(name=self.in_name))
+
+
+# =============================================================================
+def raster_copy(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
+ t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
+ nodata=None):
+ """Copy a band of raster into the output file.
+
+ Function copied from gdal_merge.py
+ """
+ if nodata is not None:
+ return raster_copy_with_nodata(s_fh, s_xoff, s_yoff, s_xsize, s_ysize,
+ s_band_n, t_fh, t_xoff, t_yoff, t_xsize,
+ t_ysize, t_band_n, nodata)
+
+ s_band = s_fh.GetRasterBand(s_band_n)
+ t_band = t_fh.GetRasterBand(t_band_n)
+
+ data = s_band.ReadRaster(s_xoff, s_yoff, s_xsize, s_ysize,
+ t_xsize, t_ysize, t_band.DataType)
+ t_band.WriteRaster(t_xoff, t_yoff, t_xsize, t_ysize, data, t_xsize,
+ t_ysize, t_band.DataType)
+
+ return 0
+
+
+def raster_copy_with_nodata(s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
+ t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
+ nodata):
+ """Copy a band of raster into the output file with nodata values.
+
+ Function copied from gdal_merge.py
+ """
+ try:
+ import numpy as Numeric
+ except ImportError:
+ import Numeric
+
+ s_band = s_fh.GetRasterBand(s_band_n)
+ t_band = t_fh.GetRasterBand(t_band_n)
+
+ data_src = s_band.ReadAsArray(s_xoff, s_yoff, s_xsize, s_ysize,
+ t_xsize, t_ysize)
+ data_dst = t_band.ReadAsArray(t_xoff, t_yoff, t_xsize, t_ysize)
+
+ nodata_test = Numeric.equal(data_src, nodata)
+ to_write = Numeric.choose(nodata_test, (data_src, data_dst))
+
+ t_band.WriteArray(to_write, t_xoff, t_yoff)
+
+ return 0
+
+
+class file_info:
+ """A class holding information about a GDAL file.
+
+ Class copied from gdal_merge.py
+
+ :param str filename: Name of file to read.
+
+ :return: 1 on success or 0 if the file can't be opened.
+ """
+
+ def init_from_name(self, filename):
+ """Initialize file_info from filename"""
+ fh = gdal.Open(filename)
+ if fh is None:
+ return 0
+
+ self.filename = filename
+ self.bands = fh.RasterCount
+ self.xsize = fh.RasterXSize
+ self.ysize = fh.RasterYSize
+ self.band_type = fh.GetRasterBand(1).DataType
+ self.block_size = fh.GetRasterBand(1).GetBlockSize()
+ self.projection = fh.GetProjection()
+ self.geotransform = fh.GetGeoTransform()
+ self.ulx = self.geotransform[0]
+ self.uly = self.geotransform[3]
+ self.lrx = self.ulx + self.geotransform[1] * self.xsize
+ self.lry = self.uly + self.geotransform[5] * self.ysize
+
+ meta = fh.GetMetadata()
+ if '_FillValue' in meta.keys():
+ self.fill_value = meta['_FillValue']
+ else:
+ self.fill_value = None
+
+ ct = fh.GetRasterBand(1).GetRasterColorTable()
+ if ct is not None:
+ self.ct = ct.Clone()
+ else:
+ self.ct = None
+
+ return 1
+
+ def copy_into(self, t_fh, s_band=1, t_band=1, nodata_arg=None):
+ """Copy this files image into target file.
+
+ This method will compute the overlap area of the file_info objects
+ file, and the target gdal.Dataset object, and copy the image data
+ for the common window area. It is assumed that the files are in
+ a compatible projection. no checking or warping is done. However,
+ if the destination file is a different resolution, or different
+ image pixel type, the appropriate resampling and conversions will
+ be done (using normal GDAL promotion/demotion rules).
+
+ :param t_fh: gdal.Dataset object for the file into which some or all
+ of this file may be copied.
+ :param s_band:
+ :param t_band:
+ :param nodata_arg:
+
+ :return: 1 on success (or if nothing needs to be copied), and zero one
+ failure.
+
+ """
+ t_geotransform = t_fh.GetGeoTransform()
+ t_ulx = t_geotransform[0]
+ t_uly = t_geotransform[3]
+ t_lrx = t_geotransform[0] + t_fh.RasterXSize * t_geotransform[1]
+ t_lry = t_geotransform[3] + t_fh.RasterYSize * t_geotransform[5]
+
+ # figure out intersection region
+ tgw_ulx = max(t_ulx, self.ulx)
+ tgw_lrx = min(t_lrx, self.lrx)
+ if t_geotransform[5] < 0:
+ tgw_uly = min(t_uly, self.uly)
+ tgw_lry = max(t_lry, self.lry)
+ else:
+ tgw_uly = max(t_uly, self.uly)
+ tgw_lry = min(t_lry, self.lry)
+
+ # do they even intersect?
+ if tgw_ulx >= tgw_lrx:
+ return 1
+ if t_geotransform[5] < 0 and tgw_uly <= tgw_lry:
+ return 1
+ if t_geotransform[5] > 0 and tgw_uly >= tgw_lry:
+ return 1
+
+ # compute target window in pixel coordinates.
+ tw_xoff = int((tgw_ulx - t_geotransform[0]) / t_geotransform[1] + 0.1)
+ tw_yoff = int((tgw_uly - t_geotransform[3]) / t_geotransform[5] + 0.1)
+ tw_xsize = int((tgw_lrx - t_geotransform[0]) / t_geotransform[1] + 0.5) - tw_xoff
+ tw_ysize = int((tgw_lry - t_geotransform[3]) / t_geotransform[5] + 0.5) - tw_yoff
+
+ if tw_xsize < 1 or tw_ysize < 1:
+ return 1
+
+ # Compute source window in pixel coordinates.
+ sw_xoff = int((tgw_ulx - self.geotransform[0]) / self.geotransform[1])
+ sw_yoff = int((tgw_uly - self.geotransform[3]) / self.geotransform[5])
+ sw_xsize = int((tgw_lrx - self.geotransform[0])
+ / self.geotransform[1] + 0.5) - sw_xoff
+ sw_ysize = int((tgw_lry - self.geotransform[3])
+ / self.geotransform[5] + 0.5) - sw_yoff
+
+ if sw_xsize < 1 or sw_ysize < 1:
+ return 1
+
+ # Open the source file, and copy the selected region.
+ s_fh = gdal.Open(self.filename)
+
+ return \
+ raster_copy(s_fh, sw_xoff, sw_yoff, sw_xsize, sw_ysize, s_band,
+ t_fh, tw_xoff, tw_yoff, tw_xsize, tw_ysize, t_band,
+ nodata_arg)
+
+
+class createMosaicGDAL:
+ """A class to mosaic modis data from hdf to GDAL formats using GDAL
+
+ :param list hdfnames: a list containing the name of tile to mosaic
+ :param str subset: the subset of layer to consider
+ :param str outformat: the output format to use, this parameter is
+ not used for the VRT output, supported values
+ are HDF4Image, GTiff, HFA, and maybe something
+ else not tested.
+ """
+ def __init__(self, hdfnames, subset, outformat="HDF4Image"):
+ """Function for the initialize the object"""
+ # Open source dataset
+ self.in_names = hdfnames
+ # #TODO use resolution into mosaic.
+ # self.resolution = res
+ if not subset:
+ self.subset = None
+ elif type(subset) == ListType:
+ self.subset = subset
+ elif type(subset) == StringType:
+ self.subset = subset.replace('(', '').replace(')',
+ '').strip().split()
+ else:
+ raise Exception('Type for subset parameter not supported')
+ self.driver = gdal.GetDriverByName(outformat)
+ if self.driver is None:
+ raise Exception('Format driver %s not found, pick a supported '
+ 'driver.' % outformat)
+ driverMD = self.driver.GetMetadata()
+ if 'DCAP_CREATE' not in driverMD:
+ raise Exception('Format driver %s does not support creation and'
+ ' piecewise writing.\nPlease select a format that'
+ ' does, such as GTiff (the default) or HFA (Erdas'
+ ' Imagine).' % format)
+ self._initLayers()
+ self._getUsedLayers()
+ self._names_to_fileinfos()
+
+ def _initLayers(self):
+ """Set up the variable self.layers as dictionary for each choosen
+ subset"""
+ if type(self.in_names) == ListType:
+ src_ds = gdal.Open(self.in_names[0])
+ else:
+ raise Exception("The input value should be a list of HDF files")
+ layers = src_ds.GetSubDatasets()
+ self.layers = {}
+ n = 0
+ if not self.subset:
+ self.subset = [1 for i in range(len(layers))]
+ for i in self.subset:
+ if str(i) == '1':
+ name = layers[n][0].split(':')[-1]
+ self.layers[name] = []
+ n = n + 1
+
+ def _getUsedLayers(self):
+ """Add each subset to the correct list for each input layers"""
+ for name in self.in_names:
+ src_ds = gdal.Open(name)
+ layers = src_ds.GetSubDatasets()
+ n = 0
+ for i in self.subset:
+ if str(i) == '1':
+ name = layers[n][0].split(':')[-1]
+ self.layers[name].append(layers[n][0])
+ n = n + 1
+
+ def _names_to_fileinfos(self):
+ """Translate a list of GDAL filenames, into file_info objects.
+ Returns a list of file_info objects. There may be less file_info
+ objects than names if some of the names could not be opened as GDAL
+ files.
+ """
+ self.file_infos = {}
+ for k, v in self.layers.iteritems():
+ self.file_infos[k] = []
+ for name in v:
+ fi = file_info()
+ if fi.init_from_name(name) == 1:
+ self.file_infos[k].append(fi)
+ self.file_infos
+
+ def _calculateNewSize(self):
+ """Return the new size of output raster
+
+ :return: X size, Y size and geotransform parameters
+ """
+ values = self.file_infos.values()
+ l1 = values[0][0]
+ ulx = l1.ulx
+ uly = l1.uly
+ lrx = l1.lrx
+ lry = l1.lry
+ for fi in self.file_infos[self.file_infos.keys()[0]]:
+ ulx = min(ulx, fi.ulx)
+ uly = max(uly, fi.uly)
+ lrx = max(lrx, fi.lrx)
+ lry = min(lry, fi.lry)
+ psize_x = l1.geotransform[1]
+ psize_y = l1.geotransform[5]
+
+ geotransform = [ulx, psize_x, 0, uly, 0, psize_y]
+ xsize = int((lrx - ulx) / geotransform[1] + 0.5)
+ ysize = int((lry - uly) / geotransform[5] + 0.5)
+ return xsize, ysize, geotransform
+
+ def write_mosaic_xml(self, prefix):
+ """Write the XML metadata file for MODIS mosaic
+
+ :param str prefix: the prefix for the XML file containing metadata
+ """
+ from parsemodis import parseModisMulti
+ import os
+ listHDF = []
+ for i in self.in_names:
+ listHDF.append(os.path.realpath(i.strip()))
+ pmm = parseModisMulti(listHDF)
+ pmm.writexml("%s.xml" % prefix)
+
+ def run(self, output):
+ """Create the mosaik
+
+ :param str output: the name of output file
+ """
+ values = self.file_infos.values()
+ l1 = values[0][0]
+ xsize, ysize, geotransform = self._calculateNewSize()
+ t_fh = self.driver.Create(output, xsize, ysize,
+ len(self.file_infos.keys()), l1.band_type)
+ if t_fh is None:
+ raise Exception('Not possibile to create dataset %s' % output)
+ return
+
+ t_fh.SetGeoTransform(geotransform)
+ t_fh.SetProjection(l1.projection)
+ i = 1
+ for names in self.file_infos.values():
+ fill = None
+ if names[0].fill_value:
+ fill = float(names[0].fill_value)
+ t_fh.GetRasterBand(i).SetNoDataValue(fill)
+ t_fh.GetRasterBand(i).Fill(fill)
+ for n in names:
+ n.copy_into(t_fh, 1, i, fill)
+ i = i + 1
+ self.write_mosaic_xml(output)
+ t_fh = None
+
+ def _calculateOffset(self, fileinfo, geotransform):
+ """Return the offset between main origin and the origin of current
+ file
+
+ :param fileinfo: a file_info object
+ :param geotransform: the geotransform parameters to keep x and y origin
+ """
+ x = abs(int((geotransform[0] - fileinfo.ulx) / geotransform[1]))
+ y = abs(int((geotransform[3] - fileinfo.uly) / geotransform[5]))
+ return x, y
+
+ def write_vrt(self, output, separate=True):
+ """Write VRT file
+
+ :param str output: the prefix of output file
+ :param bool separate: True to write a VRT file for each band, False to
+ write an unique file
+ """
+
+ def write_complex(outfile, f, geot):
+ """Write a complex source to VRT file"""
+ out.write('\t\t<ComplexSource>\n')
+ out.write('\t\t\t<SourceFilename relativeToVRT="0">{name}'
+ '</SourceFilename>\n'.format(name=f.filename.replace('"', '')))
+ out.write('\t\t\t<SourceBand>1</SourceBand>\n')
+ out.write('\t\t\t<SourceProperties RasterXSize="{x}" '
+ 'RasterYSize="{y}" DataType="{typ}" '
+ 'BlockXSize="{bx}" BlockYSize="{by}" />'
+ '\n'.format(x=f.xsize, y=f.ysize,
+ typ=f.band_type, bx=f.block_size[0],
+ by=f.block_size[1]))
+ out.write('\t\t\t<SrcRect xOff="0" yOff="0" xSize="{x}" '
+ 'ySize="{y}" />\n'.format(x=f.xsize, y=f.ysize))
+ xoff, yoff = self._calculateOffset(f, geot)
+ out.write('\t\t\t<DstRect xOff="{xoff}" yOff="{yoff}" '
+ 'xSize="{x}" ySize="{y}" />'
+ '\n'.format(xoff=xoff, yoff=yoff, x=f.xsize,
+ y=f.ysize))
+ if l1.fill_value:
+ out.write('\t\t\t<NODATA>{va}</NODATA>'
+ '\n'.format(va=f.fill_value))
+ out.write('\t\t</ComplexSource>\n')
+
+ xsize, ysize, geot = self._calculateNewSize()
+ if separate:
+ for k in self.file_infos.keys():
+ l1 = self.file_infos[k][0]
+ out = open("{pref}_{band}.vrt".format(pref=output, band=k),
+ 'w')
+ out.write('<VRTDataset rasterXSize="{x}" rasterYSize="{y}">'
+ '\n'.format(x=xsize, y=ysize))
+ out.write('\t<SRS>{proj}</SRS>\n'.format(proj=l1.projection))
+ out.write('\t<GeoTransform>{geo0}, {geo1}, {geo2}, {geo3},'
+ ' {geo4}, {geo5}</GeoTransform>'
+ '\n'.format(geo0=geot[0], geo1=geot[1], geo2=geot[2],
+ geo3=geot[3], geo4=geot[4],
+ geo5=geot[5]))
+ gtype = gdal.GetDataTypeName(l1.band_type)
+ out.write('\t<VRTRasterBand dataType="{typ}" band="1"'
+ '>\n'.format(typ=gtype))
+ if l1.fill_value:
+ out.write('\t\t<NoDataValue>{val}</NoDataValue'
+ '>\n'.format(val=l1.fill_value))
+ out.write('<ColorInterp>Gray</ColorInterp>\n')
+ for f in self.file_infos[k]:
+ write_complex(out, f, geot)
+ out.write('\t</VRTRasterBand>\n')
+ out.write('</VRTDataset>\n')
+ out.close()
+ else:
+ values = self.file_infos.values()
+ l1 = values[0][0]
+ band = 1 # the number of band
+ out = open("{pref}.vrt".format(pref=output), 'w')
+ out.write('<VRTDataset rasterXSize="{x}" rasterYSize="{y}">'
+ '\n'.format(x=xsize, y=ysize))
+ out.write('\t<SRS>{proj}</SRS>\n'.format(proj=l1.projection))
+ out.write('\t<GeoTransform>{geo0}, {geo1}, {geo2}, {geo3},'
+ ' {geo4}, {geo5}</GeoTransform>\n'.format(geo0=geot[0],
+ geo1=geot[1], geo2=geot[2], geo3=geot[3], geo4=geot[4],
+ geo5=geot[5]))
+ for k in self.file_infos.keys():
+ l1 = self.file_infos[k][0]
+ out.write('\t<VRTRasterBand dataType="{typ}" band="{n}"'
+ '>\n'.format(typ=gdal.GetDataTypeName(l1.band_type),
+ n=band))
+ if l1.fill_value:
+ out.write('\t\t<NoDataValue>{val}</NoDataValue>\n'.format(
+ val=l1.fill_value))
+ out.write('\t\t<ColorInterp>Gray</ColorInterp>\n')
+ for f in self.file_infos[k]:
+ write_complex(out, f, geot)
+ out.write('\t</VRTRasterBand>\n')
+ band += 1
+ out.write('</VRTDataset>\n')
+ out.close()
Modified: grass-addons/grass7/raster/r.modis/libmodis/downmodis.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/downmodis.py 2015-10-19 01:03:34 UTC (rev 66525)
+++ grass-addons/grass7/raster/r.modis/libmodis/downmodis.py 2015-10-19 13:43:18 UTC (rev 66526)
@@ -2,8 +2,11 @@
# class to download modis data
#
# (c) Copyright Luca Delucchi 2010
+# (c) Copyright Logan C Byers 2014
# Authors: Luca Delucchi
+# Logan C Byers
# Email: luca dot delucchi at iasma dot it
+# loganbyers at ku.edu
#
##################################################################
#
@@ -18,8 +21,25 @@
# See the GNU General Public License for more details.
#
##################################################################
+"""Module to download MODIS HDF files from NASA repository.
+It supports both FTP and HTTP repositories
-from datetime import date, timedelta
+Classes:
+
+* :class:`modisHtmlParser`
+* :class:`downModis`
+
+Functions:
+
+* :func:`urljoin`
+* :func:`getNewerVersion`
+* :func:`str2date`
+
+"""
+from __future__ import print_function
+
+from datetime import date
+from datetime import timedelta
import os
import glob
import logging
@@ -28,64 +48,82 @@
import ftplib
try:
import requests
-except:
- pass
-
-import urllib2
+except ImportError:
+ import urllib2
from HTMLParser import HTMLParser
import re
+global GDAL
+
try:
import osgeo.gdal as gdal
+ GDAL = True
except ImportError:
try:
import gdal
+ GDAL = True
except ImportError:
- raise('Python GDAL library not found, please install python-gdal')
+ GDAL = False
+ print('WARNING: Python GDAL library not found, please install it to'
+ ' check data downloaded with pyModis')
+# setup gdal
+if GDAL:
+ gdal.UseExceptions()
+ gdalDriver = gdal.GetDriverByName('HDF4')
+ if not gdalDriver:
+ GDAL = False
+ print("GDAL installation has no support for HDF4, please update GDAL")
def urljoin(*args):
- """
- Joins given arguments into a url. Trailing but not leading slashes are
+ """Joins given arguments into a url. Trailing but not leading slashes are
stripped for each argument.
http://stackoverflow.com/a/11326230
+
+ :return: a string
"""
return "/".join(map(lambda x: str(x).rstrip('/'), args))
def getNewerVersion(oldFile, newFile):
- """ Return newer version of a file
+ """Check two files to determine which is newer
- oldFile = one of the two similar file
+ :param str oldFile: one of the two similar files
+ :param str newFile: one of the two similar files
- newFile = one of the two similar file
+ :return: the name of newer file
"""
- oldFileSplit = oldFile.split('.')
- newFileSplit = newFile.split('.')
- if oldFileSplit[4] > newFileSplit[4]:
+ # get the processing date (YYYYDDDHHMMSS) from the file strings
+ if oldFile.split('.')[4] > newFile.split('.')[4]:
return oldFile
else:
return newFile
-def str2date(strin):
- """Return a date object from a string
+def str2date(datestring):
+ """Convert to datetime.date object from a string
- string = text string to return date (2012-10-04)
+ :param str datestring string with format (YYYY-MM-DD)
+ :return: a datetime.date object representing datestring
"""
- todaySplit = strin.split('-')
- return date(int(todaySplit[0]), int(todaySplit[1]), int(todaySplit[2]))
+ if '-' in datestring:
+ stringSplit = datestring.split('-')
+ elif '.' in datestring:
+ stringSplit = datestring.split('.')
+ elif ' ' in datestring:
+ stringSplit = datestring.split(' ')
+ return date(int(stringSplit[0]), int(stringSplit[1]), int(stringSplit[2]))
class modisHtmlParser(HTMLParser):
- """A class to parse HTML"""
+ """A class to parse HTML
+
+ :param fh: content of http request
+ """
def __init__(self, fh):
- """
- {fh} must be astring returned by requests.content or
- urllib2.urlopen().read()
- """
+ """Function to initialize the object"""
HTMLParser.__init__(self)
self.fileids = []
self.feed(fh)
@@ -96,96 +134,92 @@
self.fileids.append(attrD['href'].replace('/', ''))
def get_all(self):
- """ Return everything """
+ """Return everything"""
return self.fileids
def get_dates(self):
- """ Return a list of directories with date """
+ """Return a list of directories with date"""
regex = re.compile('(\d{4})[/.-](\d{2})[/.-](\d{2})$')
return [elem for elem in self.fileids if regex.match(elem)]
def get_tiles(self, prod, tiles, jpeg=False):
- """ Return a list of file to download """
+ """Return a list of files to download
+
+ :param str prod: the code of MODIS product that we are going to
+ analyze
+ :param list tiles: the list of tiles to consider
+ :param bool jpeg: True to also check for jpeg data
+ """
finalList = []
for i in self.fileids:
+ # distinguish jpg from hdf by where the tileID is within the string
+ # jpgs have the tileID at index 3, hdf have tileID at index 2
name = i.split('.')
- # distinguish jpeg files from hdf files by the number of index
- # where find the tile index
+ # if product is not in the filename, move to next filename in list
if not name.count(prod):
continue
+ # if tiles are not specified and the file is not a jpg, add to list
if not tiles and not (name.count('jpg') or name.count('BROWSE')):
finalList.append(i)
- #is a jpeg of tiles number
+ # if tiles are specified
if tiles:
+ # if a tileID is at index 3 and jpgs are to be downloaded
if tiles.count(name[3]) == 1 and jpeg:
finalList.append(i)
- #is a hdf of tiles number
+ # if a tileID is at in index 2, it is known to be HDF
elif tiles.count(name[2]) == 1:
finalList.append(i)
return finalList
class downModis:
- """A class to download MODIS data from NASA FTP repository"""
- def __init__(self,
- destinationFolder,
- password=None,
- user="anonymous",
- url="http://e4ftl01.cr.usgs.gov",
- tiles=None,
- path="MOLT",
- product="MOD11A1.005",
- today=None,
- enddate=None,
- delta=10,
- jpg=False,
- debug=False,
- timeout=30
- ):
- """Initialization function :
+ """A class to download MODIS data from NASA FTP or HTTP repositories
- destinationFolder = where the files will be stored
+ :param str destinationFolder: where the files will be stored
+ :param str password: the password, it should be your email address
+ to connect to a FTP server. Do not use this
+ variable if the server is an HTTP server
+ :param str user: the user name, by default 'anonymous', used to
+ connect to an FTP server. Do not use this variable
+ if the server is an HTTP server
+ :param str url: the base url from where to download the MODIS data,
+ it can be FTP or HTTP but it has to start with
+ 'ftp://' or 'http://'
+ :param str path: the directory where the data that you want to
+ download are stored on the FTP server. For HTTP
+ requests, this is the part of the url between the 'url'
+ parameter and the 'product' parameter.
+ :param str product: the code of the product to download, the code
+ should be idential to the one of the url
+ :param str tiles: a set of tiles to be downloaded, None == all tiles.
+ This can be passed as a string of tileIDs separated
+ by commas, or as a list of individual tileIDs
+ :param str today: the day to start downloading; in order to pass a
+ date different from today use the format YYYY-MM-DD
+ :param str enddate: the day to end downloading; in order to pass a
+ date use the format YYYY-MM-DD. This day must be
+ before the 'today' parameter. Downloading happens
+ in reverse order (currently)
+ :param int delta: timelag i.e. the number of days starting from
+ today backwards. Will be overwritten if
+ 'enddate' is specifed during instantiation
+ :param bool jpeg: set to True if you want to download the JPG overview
+ file in addition to the HDF
+ :param bool debug: set to True if you want to obtain debug information
+ :param int timeout: Timeout value for HTTP server (seconds)
+ :param bool checkgdal: variable to set the GDAL check
+ """
+ def __init__(self, destinationFolder, password=None, user="anonymous",
+ url="http://e4ftl01.cr.usgs.gov", tiles=None, path="MOLT",
+ product="MOD11A1.005", today=None, enddate=None, delta=10,
+ jpg=False, debug=False, timeout=30, checkgdal=True):
+ """Function to initialize the object"""
- password = the password, it should be your email address to
- connect to a FTP server.
- Do not use this variable if the server is a HTTP server
-
- user = the user name, by default 'anonymous', used to connect
- to a FTP server.
- Do not use this variable if the server is a HTTP server
-
- url = the url from where to download the MODIS data, it can be FTP or
- HTTP and it has to start with http:// or ftp://
-
- path = the directory where the data that you want to download are
- stored on the FTP server
-
- product = the code of product to download, the code should be
- idential to the one of the url
-
- tiles = a list of tiles to be downloaded, None == all tiles
-
- today = the day to start downloading; in order to pass a date
- different from today use the format YYYY-MM-DD
-
- enddate = the day to end downloading; in order to pass a date
- use the format YYYY-MM-DD
-
- delta = timelag i.e. the number of days starting from today
- backwards
-
- jpeg = set to True if you want to download also the JPG overview file
-
- debug = set to True if you want to obtain debug information
-
- timeout = Timeout value for HTTP server
- """
-
- # url modis
- if 'ftp' in url:
+ # prepare the base url and set the url type (ftp/http)
+ if 'ftp://' in url:
self.url = url.replace('ftp://', '').rstrip('/')
self.urltype = 'ftp'
- elif 'http' in url:
+ elif 'http://' in url:
self.url = url
self.urltype = 'http'
else:
@@ -194,65 +228,71 @@
self.user = user
# password for download using ftp
self.password = password
- # the product
+ # the product (product_code.004 or product_cod.005)
self.product = product
self.product_code = product.split('.')[0]
- # directory where data are collected
+ # url directory where data are located
self.path = urljoin(path, self.product)
# tiles to downloads
- if tiles:
+ if type(tiles) == str:
self.tiles = tiles.split(',')
- else:
+ else: # tiles are list, tuple, or None
self.tiles = tiles
# set destination folder
if os.access(destinationFolder, os.W_OK):
self.writeFilePath = destinationFolder
else:
- raise IOError("Folder to store downloaded files does not exist" \
- + " or is not writeable")
+ raise Exception("Folder to store downloaded files does not exist"
+ " or is not writeable")
# return the name of product
if len(self.path.split('/')) == 2:
self.product = self.path.split('/')[1]
elif len(self.path.split('/')) == 3:
self.product = self.path.split('/')[2]
- # write a file with the name of file downloaded
- self.filelist = open(os.path.join(self.writeFilePath, 'listfile' \
- + self.product + '.txt'), 'w')
- # set jpg download
+ # write a file with the name of file to be downloaded
+ self.filelist = open(os.path.join(self.writeFilePath,
+ 'listfile{pro}.txt'.format(pro=self.product)),
+ 'w')
+ # set if to download jpgs
self.jpeg = jpg
- # today
+ # today, or the last day in the download series chronologically
self.today = today
- # force the last day
+ # chronologically the first day in the download series
self.enday = enddate
- # delta of days
+ # default number of days to consider if enddate not specified
self.delta = delta
# status of tile download
self.status = True
# for debug, you can download only xml files
self.debug = debug
# for logging
- log_filename = os.path.join(self.writeFilePath, 'modis' \
- + self.product + '.log')
+ log_filename = os.path.join(self.writeFilePath,
+ 'modis{pro}.log'.format(pro=self.product))
log_format = '%(asctime)s - %(levelname)s - %(message)s'
- logging.basicConfig(filename=log_filename, level=logging.DEBUG, \
- format=log_format)
+ logging.basicConfig(filename=log_filename, level=logging.DEBUG,
+ format=log_format)
+ # global connection attempt counter
self.nconnection = 0
+ # timeout for HTTP connection before failing (seconds)
self.timeout = timeout
+ # files within the directory where data will be saved
self.fileInPath = []
- # add all files in the directory where we will save new modis data
for f in os.listdir(self.writeFilePath):
if os.path.isfile(os.path.join(self.writeFilePath, f)):
self.fileInPath.append(f)
- gdal.UseExceptions()
- gdalDriver = gdal.GetDriverByName('HDF4')
- if not gdalDriver:
- raise IOError("GDAL installation has no support for HDF4, please update GDAL")
+ global GDAL
+ if not GDAL and checkgdal:
+ logging.warning("WARNING: Python GDAL library is not found")
+ elif GDAL and not checkgdal:
+ GDAL = False
def removeEmptyFiles(self):
- """Check if some file has size equal to 0"""
+ """Function to remove files in the download directory that have
+ filesize equal to 0
+ """
year = str(date.today().year)
- pref = self.product.split('.')[0]
- files = glob.glob1(self.writeFilePath, '%s.A%s*' % (pref, year))
+ prefix = self.product.split('.')[0]
+ files = glob.glob1(self.writeFilePath, '%s.A%s*' % (prefix, year))
for f in files:
fil = os.path.join(self.writeFilePath, f)
if os.path.getsize(fil) == 0:
@@ -261,19 +301,23 @@
def connect(self, ncon=20):
"""Connect to the server and fill the dirData variable
- ncon = maximum number of attempts to connect to the HTTP server
- before failing
+ :param int ncon: maximum number of attempts to connect to the HTTP
+ server before failing
"""
if self.urltype == 'ftp':
self._connectFTP(ncon)
elif self.urltype == 'http':
self._connectHTTP(ncon)
+ if len(self.dirData) == 0:
+ raise Exception("There are some troubles with the server. "
+ "The directory seems to be empty")
def _connectHTTP(self, ncon=20):
- """ Connect to HTTP server, create a list of directories for all days
+ """Connect to HTTP server, create a list of directories for all days
- ncon = maximum number of attempts to connect to the HTTP server
- before failing
+ :param int ncon: maximum number of attempts to connect to the HTTP
+ server before failing. If ncon < 0, connection
+ attempts are unlimited in number
"""
self.nconnection += 1
try:
@@ -288,15 +332,15 @@
self.dirData.reverse()
except:
logging.error('Error in connection')
- if self.nconnection <= ncon:
+ if self.nconnection <= ncon or ncon < 0:
self._connectHTTP()
def _connectFTP(self, ncon=20):
- """ Set connection to ftp server, move to path where data are stored
- and create a list of directories for all days
+ """Set connection to ftp server, move to path where data are stored,
+ and create a list of directories for all days
- ncon = maximum number of attempts to connect to the FTP server
- before failing
+ :param int ncon: maximum number of attempts to connect to the FTP
+ server before failing.
"""
self.nconnection += 1
@@ -311,58 +355,63 @@
self.ftp.dir(self.dirData.append)
# reverse order of data for have first the nearest to today
self.dirData.reverse()
- # check if dirData contain only directory, delete all files
+ # ensure dirData contains only directories, remove all references to files
self.dirData = [elem.split()[-1] for elem in self.dirData if elem.startswith("d")]
- if self.debug == True:
- logging.debug("Open connection %s" % self.url)
+ if self.debug:
+ logging.debug("Open connection {url}".format(url=self.url))
except (EOFError, ftplib.error_perm), e:
- logging.error('Error in connection: %s' % e)
+ logging.error('Error in connection: {err}'.format(err=e))
if self.nconnection <= ncon:
self._connectFTP()
def closeFTP(self):
- """ Close ftp connection """
+ """Close ftp connection and close the file list document"""
self.ftp.quit()
self.closeFilelist()
- if self.debug == True:
- logging.debug("Close connection %s" % self.url)
+ if self.debug:
+ logging.debug("Close connection {url}".format(url=self.url))
def closeFilelist(self):
- """ Function to close the file where write the downloaded files """
+ """Function to close the file list of where the files are downloaded"""
self.filelist.close()
def setDirectoryIn(self, day):
- """ Enter in the directory of the day """
+ """Enter into the file directory of a specified day
+
+ :param str day: a string representing a day in format YYYY.MM.DD
+ """
try:
self.ftp.cwd(day)
except (ftplib.error_reply, socket.error), e:
- logging.error("Error %s entering in directory %s" % e, day)
+ logging.error("Error {err} entering in directory "
+ "{name}".format(err=e, name=day))
self.setDirectoryIn(day)
def setDirectoryOver(self):
- """ Come back to old path """
+ """Move up within the file directory"""
try:
self.ftp.cwd('..')
except (ftplib.error_reply, socket.error), e:
- logging.error("Error %s when try to come back" % e)
+ logging.error("Error {err} when trying to come back".format(err=e))
self.setDirectoryOver()
def _getToday(self):
- """Return the first day for start to download"""
- if self.today == None:
- # set today variable to today
+ """Set the dates for the start and end of downloading"""
+ if self.today is None:
+ # set today variable from datetime.date method
self.today = date.today()
- else:
- # set today variable to data pass from user
+ elif type(self.today) == str:
+ # set today variable from string data passed by user
self.today = str2date(self.today)
- # set enday variable to data
- if self.enday != None:
+ # set enday variable to data passed by user
+ if type(self.enday) == str:
self.enday = str2date(self.enday)
+ # set delta
if self.today and self.enday:
if self.today < self.enday:
- raise IOError("The first day should be newer then end date")
- D = self.today - self.enday
- self.delta = D.days
+ self.today, self.enday = self.enday, self.today
+ delta = self.today - self.enday
+ self.delta = abs(delta.days) + 1
def getListDays(self):
"""Return a list of all selected days"""
@@ -381,12 +430,13 @@
days = self.dirData[today_index:][:self.delta]
# this is useful for 8/16 days data, delta could download more images
# that you want
- if self.enday != None:
+ if self.enday is not None:
enday_s = self.enday.strftime("%Y.%m.%d")
delta = 0
+ # make a full cycle from the last index and find
# it make a for cicle from the last value and find the internal
# delta to remove file outside temporaly range
- for i in range(-(len(days)), 0):
+ for i in range(0, len(days)):
if days[i] < enday_s:
break
else:
@@ -400,11 +450,13 @@
return self.dirData
def getFilesList(self, day=None):
- """ Creates a list of files to download, it is possible to choose to
- download also the JPG overview files or only the HDF files
+ """Returns a list of files to download. HDF and XML files are
+ downloaded by default. JPG files will be downloaded if
+ self.jpeg == True.
- day = the date of data
+ :param str day: the date of data in format YYYY.MM.DD
+ :return: a list of files to download for the day
"""
if self.urltype == 'http':
return self._getFilesListHTTP(day)
@@ -412,46 +464,47 @@
return self._getFilesListFTP()
def _getFilesListHTTP(self, day):
- """ Creates a list of files to download from http server, it is
- possible to choose to download also the JPG overview files or
- only the HDF files
+ """Returns a list of files to download from http server, which will
+ be HDF and XML files, and optionally JPG files if specified by
+ self.jpeg
- day = the date of data
-
+ :param str day: the date of data in format YYYY.MM.DD
"""
- # return the file's list inside the directory of each day
+ # return the files list inside the directory of each day
try:
url = urljoin(self.url, self.path, day)
- if self.debug == True:
- logging.debug("The url is: %s" % url)
+ if self.debug:
+ logging.debug("The url is: {url}".format(url=url))
try:
http = modisHtmlParser(requests.get(url,
timeout=self.timeout).content)
except:
http = modisHtmlParser(urllib2.urlopen(url,
timeout=self.timeout).read())
- # download also jpeg
+ # download JPG files also
if self.jpeg:
- # finallist is ugual to all file with jpeg file
+ # if tiles not specified, download all files
if not self.tiles:
finalList = http.get_all()
- # finallist is ugual to tiles file with jpeg file
+ # if tiles specified, download all files with jpegs
else:
finalList = http.get_tiles(self.product_code,
self.tiles, jpeg=True)
- # not download jpeg
+ # if JPG files should not be downloaded, get only HDF and XML
else:
finalList = http.get_tiles(self.product_code, self.tiles)
- if self.debug == True:
- logging.debug("The number of file to download is: %i" % len(finalList))
+ if self.debug:
+ logging.debug("The number of file to download is: "
+ "{num}".format(num=len(finalList)))
return finalList
except (socket.error), e:
- logging.error("Error %s when try to receive list of files" % e)
+ logging.error("Error {err} when try to receive list of "
+ "files".format(err=e))
self._getFilesListHTTP(day)
def _getFilesListFTP(self):
- """ Create a list of files to download from FTP server, it is possible
- choose to download also the JPG overview files or only the HDF files
+ """Create a list of files to download from FTP server, it is possible
+ choose to download also the JPG overview files or only the HDF files
"""
def cicle_file(jpeg=False):
"""Check the type of file"""
@@ -463,11 +516,11 @@
if not self.tiles and not (name.count('jpg') or
name.count('BROWSE')):
finalList.append(i)
- #is a jpeg of tiles number
+ # is a jpeg of tiles number
if self.tiles:
if self.tiles.count(name[3]) == 1 and jpeg:
finalList.append(i)
- #is a hdf of tiles number
+ # is a hdf of tiles number
elif self.tiles.count(name[2]) == 1:
finalList.append(i)
return finalList
@@ -486,33 +539,38 @@
# not download jpeg
else:
finalList = cicle_file()
- if self.debug == True:
- logging.debug("The number of file to download is: %i" % len(finalList))
+ if self.debug:
+ logging.debug("The number of file to download is: "
+ "{num}".format(num=len(finalList)))
return finalList
except (ftplib.error_reply, socket.error), e:
- logging.error("Error %s when trying to receive list of files" % e)
+ logging.error("Error {err} when trying to receive list of "
+ "files".format(err=e))
self._getFilesListFTP()
- def checkDataExist(self, listNewFile, move=0):
- """ Check if a file already exists in the directory of download
+ def checkDataExist(self, listNewFile, move=False):
+ """Check if a file already exists in the local download directory
- listNewFile = list of all files, returned by getFilesList function
-
- move = it is useful to know if a function is called from download
- or move function
+ :param list listNewFile: list of all files, returned by getFilesList
+ function
+ :param bool move: it is useful to know if a function is called from
+ download or move function
+ :return: list of files to download
"""
# different return if this method is used from downloadsAllDay() or
# moveFile()
- if move == 0:
+ if not move:
listOfDifferent = list(set(listNewFile) - set(self.fileInPath))
- elif move == 1:
+ elif move:
listOfDifferent = list(set(self.fileInPath) - set(listNewFile))
return listOfDifferent
def checkFile(self, filHdf):
"""Check by using GDAL to be sure that the download went ok
- filHdf = name of the HDF file to check
+ :param str filHdf: name of the HDF file to check
+
+ :return: 0 if file is correct, 1 for error
"""
try:
gdal.Open(filHdf)
@@ -522,13 +580,11 @@
return 1
def downloadFile(self, filDown, filHdf, day):
- """Download the single file
+ """Download a single file
- filDown = name of the file to download
-
- filHdf = name of the file to write to
-
- day = the day in format YYYY.MM.DD
+ :param str filDown: name of the file to download
+ :param str filHdf: name of the file to write to
+ :param str day: the day in format YYYY.MM.DD
"""
if self.urltype == 'http':
self._downloadFileHTTP(filDown, filHdf, day)
@@ -536,78 +592,80 @@
self._downloadFileFTP(filDown, filHdf)
def _downloadFileHTTP(self, filDown, filHdf, day):
- """Download the single file from http server
+ """Download a single file from the http server
- filDown = name of the file to download
-
- filSave = name of the file to write to
-
- day = the day in format YYYY.MM.DD
+ :param str filDown: name of the file to download
+ :param str filHdf: name of the file to write to
+ :param str day: the day in format YYYY.MM.DD
"""
filSave = open(filHdf, "wb")
- orig_size = None
- try:
- try:
+ try: # download and write the file
+ try: # use request module
http = requests.get(urljoin(self.url, self.path, day, filDown))
orig_size = http.headers['content-length']
filSave.write(http.content)
- except:
+ except: # use urllib2 module
http = urllib2.urlopen(urljoin(self.url, self.path, day,
filDown))
orig_size = http.headers['content-length']
filSave.write(http.read())
- filSave.close()
- #if it have an error it try to download again the file
+
+ # if local file has an error, try to download the file again
except:
- logging.error("Cannot download %s. Retrying.." % filDown)
+ logging.error("Cannot download {name}. "
+ "Retrying..".format(name=filDown))
filSave.close()
os.remove(filSave.name)
self._downloadFileHTTP(filDown, filHdf, day)
+ filSave.close()
transf_size = os.path.getsize(filSave.name)
- if orig_size and int(orig_size) == int(transf_size):
+ if int(orig_size) == int(transf_size):
+ # if no xml file, delete the HDF and redownload
if filHdf.find('.xml') == -1:
- if self.checkFile(filHdf):
+ test = False
+ if GDAL:
+ test = self.checkFile(filHdf)
+ if test:
os.remove(filSave.name)
self._downloadFileHTTP(filDown, filHdf, day)
else:
- self.filelist.write("%s\n" % filDown)
- if self.debug == True:
- logging.debug("File %s downloaded correctly" % filDown)
+ self.filelist.write("{name}\n".format(name=filDown))
+ if self.debug:
+ logging.debug("File {name} downloaded "
+ "correctly".format(name=filDown))
return 0
- else:
- self.filelist.write("%s\n" % filDown)
- if self.debug == True:
- logging.debug("File %s downloaded correctly" % filDown)
+ else: # xml exists
+ self.filelist.write("{name}\n".format(name=filDown))
+ if self.debug:
+ logging.debug("File {name} downloaded "
+ "correctly".format(name=filDown))
return 0
+ # if filesizes are different, delete and try again
else:
- if not orig_size:
- logging.warning("Different size for file %s - original data:"\
- " None, downloaded: %s" % (filDown,
- transf_size))
- else:
- logging.warning("Different size for file %s - original data:"\
- " %s, downloaded: %s" % (filDown, orig_size,
- transf_size))
+ logging.warning("Different size for file {name} - original data: "
+ "{orig}, downloaded: {down}".format(name=filDown,
+ orig=orig_size,
+ down=transf_size))
os.remove(filSave.name)
self._downloadFileHTTP(filDown, filHdf, day)
def _downloadFileFTP(self, filDown, filHdf):
- """Download the single file from ftp server
+ """Download a single file from ftp server
- filDown = name of the file to download
-
- filSave = name of the file to write to
+ :param str filDown: name of the file to download
+ :param str filHdf: name of the file to write to
"""
filSave = open(filHdf, "wb")
- try:
+ try: # transfer file from ftp
self.ftp.retrbinary("RETR " + filDown, filSave.write)
- self.filelist.write("%s\n" % filDown)
- if self.debug == True:
- logging.debug("File %s downloaded" % filDown)
- #if it have an error it try to download again the file
- except (ftplib.error_reply, socket.error, ftplib.error_temp, EOFError), e:
- logging.error("Cannot download %s, the error was '%s'. Retrying..." % (
- filDown, e))
+ self.filelist.write("{name}\n".format(name=filDown))
+ if self.debug:
+ logging.debug("File {name} downloaded".format(name=filDown))
+ # if error during download process, try to redownload the file
+ except (ftplib.error_reply, socket.error, ftplib.error_temp,
+ EOFError), e:
+ logging.error("Cannot download {name}, the error was '{err}'. "
+ "Retrying...".format(name=filDown, err=e))
filSave.close()
os.remove(filSave.name)
try:
@@ -621,98 +679,116 @@
if orig_size == transf_size:
return 0
else:
- logging.warning("Different size for file %s - original data: %s," \
- " downloaded: %s" % (filDown, orig_size,
- transf_size))
+ logging.warning("Different size for file {name} - original data: "
+ "{orig}, downloaded: {down}".format(name=filDown,
+ orig=orig_size,
+ down=transf_size))
os.remove(filSave.name)
self._downloadFileFTP(filDown, filHdf)
def dayDownload(self, day, listFilesDown):
- """ Downloads tiles are in files_hdf_consider
+ """Downloads tiles for the selected day
- listFilesDown = list of the files to download, returned by
- checkDataExist function
+ :param str day: the day in format YYYY.MM.DD
+ :param list listFilesDown: list of the files to download, returned
+ by checkDataExist function
"""
# for each file in files' list
for i in listFilesDown:
fileSplit = i.split('.')
- filePrefix = "%s.%s.%s.%s" % (fileSplit[0], fileSplit[1],
- fileSplit[2], fileSplit[3])
+ filePrefix = "{a}.{b}.{c}.{d}".format(a=fileSplit[0],
+ b=fileSplit[1],
+ c=fileSplit[2],
+ d=fileSplit[3])
- # check data exists in the return directory
- oldFile = glob.glob1(self.writeFilePath, filePrefix + "*" \
- + fileSplit[-1])
+ # check if this file already exists in the save directory
+ oldFile = glob.glob1(self.writeFilePath, filePrefix + "*"
+ + fileSplit[-1])
numFiles = len(oldFile)
+ # if it doesn't exist
if numFiles == 0:
file_hdf = os.path.join(self.writeFilePath, i)
+ # if one does exist
elif numFiles == 1:
- # check the version of file
+ # check the version of file, delete local file if it is older
fileDown = getNewerVersion(oldFile[0], i)
if fileDown != oldFile[0]:
os.remove(os.path.join(self.writeFilePath, oldFile[0]))
file_hdf = os.path.join(self.writeFilePath, fileDown)
elif numFiles > 1:
- logging.error("There are to many files for %s" % i)
+ logging.error("There are to many files for "
+ "{name}".format(name=i))
if numFiles == 0 or (numFiles == 1 and fileDown != oldFile[0]):
self.downloadFile(i, file_hdf, day)
def downloadsAllDay(self, clean=False, allDays=False):
- """Download the single file
+ """Download all requested days
- filDown = name of the file to download
-
- filSave = name of the file to write to
+ :param bool clean: if True remove the empty files, they could have
+ some problems in the previous download
+ :param bool allDays: download all passable days
"""
- #return the days to download
if clean:
self.removeEmptyFiles()
+ # get the days to download
if allDays:
days = self.getAllDays()
else:
days = self.getListDays()
- if self.debug == True:
- logging.debug("The number of days to download is: %i" % len(days))
+ # log the days to download
+ if self.debug:
+ logging.debug("The number of days to download is: "
+ "{num}".format(num=len(days)))
+ # download the data
if self.urltype == 'http':
- self._downloadsAllDayHTTP(days)
+ self._downloadAllDaysHTTP(days)
elif self.urltype == 'ftp':
- self._downloadAllDayFTP(days)
+ self._downloadAllDaysFTP(days)
- def _downloadsAllDayHTTP(self, days):
- """ Downloads all the tiles considered from HTTP server"""
+ def _downloadAllDaysHTTP(self, days):
+ """Downloads all the tiles considered from HTTP server
- #for each day
+ :param list days: the list of days to download
+ """
+ # for each day
for day in days:
- #obtain list of all files
+ # obtain list of all files
listAllFiles = self.getFilesList(day)
- #obtain list of files to download
+ # filter files based on local files in save directory
listFilesDown = self.checkDataExist(listAllFiles)
- #download files for a day
+ # download files for a day
self.dayDownload(day, listFilesDown)
self.closeFilelist()
- if self.debug == True:
+ if self.debug:
logging.debug("Download terminated")
return 0
- def _downloadsAllDayFTP(self, days):
- """ Downloads all the tiles considered from FTP server"""
- #for each day
+ def _downloadAllDaysFTP(self, days):
+ """Downloads all the tiles considered from FTP server
+
+ :param list days: the list of days to download
+ """
+ # for each day
for day in days:
- #enter in the directory of day
+ # enter in the directory of day
self.setDirectoryIn(day)
- #obtain list of all files
+ # obtain list of all files
listAllFiles = self.getFilesList()
- #obtain list of files to download
+ # filter files based on local files in save directory
listFilesDown = self.checkDataExist(listAllFiles)
- #download files for a day
+ # download files for a day
self.dayDownload(day, listFilesDown)
self.setDirectoryOver()
self.closeFTP()
- if self.debug == True:
+ if self.debug:
logging.debug("Download terminated")
return 0
def debugLog(self):
- """Function to create the debug file"""
+ """Function to create the debug file
+
+ :return: a Logger object to use to write debug info
+ """
# create logger
logger = logging.getLogger("PythonLibModis debug")
logger.setLevel(logging.DEBUG)
@@ -720,11 +796,11 @@
ch = logging.StreamHandler()
ch.setLevel(logging.DEBUG)
# create formatter
- formatter = logging.Formatter("%(asctime)s - %(name)s - " \
- + "%(levelname)s - %(message)s")
- # add formatter to ch
+ formatter = logging.Formatter("%(asctime)s - %(name)s - "
+ "%(levelname)s - %(message)s")
+ # add formatter to console handler
ch.setFormatter(formatter)
- # add ch to logger
+ # add console handler to logger
logger.addHandler(ch)
return logger
@@ -732,24 +808,25 @@
"""This function is useful to debug the number of days"""
logger = self.debugLog()
days = self.getListDays()
- # if lenght of list of days and the delta of day they are different
+ # if length of list of days and the delta of days are different
if len(days) != self.delta:
# for each day
for i in range(1, self.delta + 1):
- # calculate the current day
+ # calculate the current day using datetime.timedelta
delta = timedelta(days=i)
day = self.today - delta
day = day.strftime("%Y.%m.%d")
# check if day is in the days list
if day not in days:
- logger.critical("This day %s is not present on list" % day)
- # the lenght of list of days and delta are ugual
+ logger.critical("This day {day} is not present on "
+ "list".format(day=day))
+ # the length of list of days and delta are equal
else:
- logger.info("All right!!")
+ logger.info("debugDays() : getListDays() and self.delta are same "
+ "length")
def debugMaps(self):
- """This function is useful to debug the number of maps to download for
- each day"""
+ """ Prints the files to download to the debug stream"""
logger = self.debugLog()
days = self.getListDays()
for day in days:
Modified: grass-addons/grass7/raster/r.modis/libmodis/parsemodis.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/parsemodis.py 2015-10-19 01:03:34 UTC (rev 66525)
+++ grass-addons/grass7/raster/r.modis/libmodis/parsemodis.py 2015-10-19 13:43:18 UTC (rev 66526)
@@ -18,15 +18,24 @@
# See the GNU General Public License for more details.
#
##################################################################
+"""Simple class to parse MODIS metadata file, it can also write the XML
+metadata file for a mosaic.
-from datetime import *
+Classes:
+
+* :class:`parseModis`
+* :class:`parseModisMulti`
+
+"""
+
+from __future__ import print_function
import string
import os
-## lists of parameters accepted by resample MRT software
+# lists of parameters accepted by resample MRT software
# projections
-PROJ_LIST = ['AEA','GEO', 'HAM', 'IGH', 'ISIN', 'LA', 'LCC', 'MOL', 'PS',
- 'SIN','TM', 'UTM', 'MERCAT']
+PROJ_LIST = ['AEA', 'GEO', 'HAM', 'IGH', 'ISIN', 'LA', 'LCC', 'MOL', 'PS',
+ 'SIN', 'TM', 'UTM', 'MERCAT']
# resampling
RESAM_LIST = ['NEAREST_NEIGHBOR', 'BICUBIC', 'CUBIC_CONVOLUTION', 'NONE']
RESAM_LIST_SWATH = ['NN', 'BI', 'CC']
@@ -38,820 +47,908 @@
class parseModis:
- """Class to parse MODIS xml files, it can also create the parameter
- configuration file for resampling MODIS DATA with the MRT software or
- convertmodis Module
- """
+ """Class to parse MODIS xml files, it can also create the parameter
+ configuration file for resampling MODIS DATA with the MRT software or
+ convertmodis Module
- def __init__(self, filename):
- """Initialization function :
-
- filename = the name of MODIS hdf file
+ :param str filename: the name of MODIS hdf file
"""
- from xml.etree import ElementTree
- if os.path.exists(filename):
- # hdf name
- self.hdfname = filename
- else:
- raise IOError('%s does not exist' % filename)
- if os.path.exists(self.hdfname + '.xml'):
- # xml hdf name
- self.xmlname = self.hdfname + '.xml'
- else:
- raise IOError('%s does not exist' % self.hdfname + '.xml')
+ def __init__(self, filename):
+ """Function to initialize the object"""
+ from xml.etree import ElementTree
+ if os.path.exists(filename):
+ # hdf name
+ self.hdfname = filename
+ else:
+ raise Exception('{name} does not exist'.format(name=filename))
- # tif name for the output file for resample MRT software
- self.tifname = self.hdfname.replace('.hdf', '.tif')
- with open(self.xmlname) as f:
- self.tree = ElementTree.parse(f)
- # return the code of tile for conf file
- self.code = os.path.split(self.hdfname)[1].split('.')[-2]
- self.path = os.path.split(self.hdfname)[0]
+ if os.path.exists(self.hdfname + '.xml'):
+ # xml hdf name
+ self.xmlname = self.hdfname + '.xml'
+ else:
+ raise Exception('{name}.xml does not exist'.format(name=self.hdfname))
- def __str__(self):
- """Print the file without xml tags"""
- retString = ""
- try:
- for node in self.tree.iter():
- if node.text.strip() != '':
- retString = "%s = %s\n" % (node.tag, node.text)
- except:
- for node in self.tree.getiterator():
- if node.text.strip() != '':
- retString = "%s = %s\n" % (node.tag, node.text)
- return retString
+ # tif name for the output file for resample MRT software
+ self.tifname = self.hdfname.replace('.hdf', '.tif')
+ with open(self.xmlname) as f:
+ self.tree = ElementTree.parse(f)
+ # return the code of tile for conf file
+ self.code = os.path.split(self.hdfname)[1].split('.')[-2]
+ self.path = os.path.split(self.hdfname)[0]
- def getRoot(self):
- """Set the root element"""
- self.rootree = self.tree.getroot()
+ def __str__(self):
+ """Print the file without xml tags"""
+ retString = ""
+ try:
+ for node in self.tree.iter():
+ if node.text.strip() != '':
+ retString = "{tag} = {val}\n".format(tag=node.tag,
+ val=node.text)
+ except:
+ for node in self.tree.getiterator():
+ if node.text.strip() != '':
+ retString = "{tag} = {val}\n".format(tag=node.tag,
+ val=node.text)
+ return retString
- def retDTD(self):
- """Return the DTDVersion element"""
- self.getRoot()
- return self.rootree.find('DTDVersion').text
+ def getRoot(self):
+ """Set the root element"""
+ self.rootree = self.tree.getroot()
- def retDataCenter(self):
- """Return the DataCenterId element"""
- self.getRoot()
- return self.rootree.find('DataCenterId').text
+ def retDTD(self):
+ """Return the DTDVersion element"""
+ self.getRoot()
+ return self.rootree.find('DTDVersion').text
- def getGranule(self):
- """Set the GranuleURMetaData element"""
- self.getRoot()
- self.granule = self.rootree.find('GranuleURMetaData')
+ def retDataCenter(self):
+ """Return the DataCenterId element"""
+ self.getRoot()
+ return self.rootree.find('DataCenterId').text
- def retGranuleUR(self):
- """Return the GranuleUR element"""
- self.getGranule()
- return self.granule.find('GranuleUR').text
+ def getGranule(self):
+ """Set the GranuleURMetaData element"""
+ self.getRoot()
+ self.granule = self.rootree.find('GranuleURMetaData')
- def retDbID(self):
- """Return the DbID element"""
- self.getGranule()
- return self.granule.find('DbID').text
+ def retGranuleUR(self):
+ """Return the GranuleUR element"""
+ self.getGranule()
+ return self.granule.find('GranuleUR').text
- def retInsertTime(self):
- """Return the InsertTime element"""
- self.getGranule()
- return self.granule.find('InsertTime').text
+ def retDbID(self):
+ """Return the DbID element"""
+ self.getGranule()
+ return self.granule.find('DbID').text
- def retLastUpdate(self):
- """Return the LastUpdate element"""
- self.getGranule()
- return self.granule.find('LastUpdate').text
+ def retInsertTime(self):
+ """Return the InsertTime element"""
+ self.getGranule()
+ return self.granule.find('InsertTime').text
- def retCollectionMetaData(self):
- """Return the CollectionMetaData element as dictionary"""
- self.getGranule()
- collect = {}
- for i in self.granule.find('CollectionMetaData').getiterator():
- if i.text.strip() != '':
- collect[i.tag] = i.text
- return collect
+ def retLastUpdate(self):
+ """Return the LastUpdate element"""
+ self.getGranule()
+ return self.granule.find('LastUpdate').text
- def retDataFiles(self):
- """Return the DataFiles element as dictionary"""
- self.getGranule()
- collect = {}
- datafiles = self.granule.find('DataFiles')
- for i in datafiles.find('DataFileContainer').getiterator():
- if i.text.strip() != '':
- collect[i.tag] = i.text
- return collect
+ def retCollectionMetaData(self):
+ """Return the CollectionMetaData element as dictionary"""
+ self.getGranule()
+ collect = {}
+ for i in self.granule.find('CollectionMetaData').getiterator():
+ if i.text.strip() != '':
+ collect[i.tag] = i.text
+ return collect
- def retDataGranule(self):
- """Return the ECSDataGranule elements as dictionary"""
- self.getGranule()
- datagran = {}
- for i in self.granule.find('ECSDataGranule').getiterator():
- if i.text.strip() != '':
- datagran[i.tag] = i.text
- return datagran
+ def retDataFiles(self):
+ """Return the DataFiles element as dictionary"""
+ self.getGranule()
+ collect = {}
+ datafiles = self.granule.find('DataFiles')
+ for i in datafiles.find('DataFileContainer').getiterator():
+ if i.text.strip() != '':
+ collect[i.tag] = i.text
+ return collect
- def retPGEVersion(self):
- """Return the PGEVersion element"""
- self.getGranule()
- return self.granule.find('PGEVersionClass').find('PGEVersion').text
+ def retDataGranule(self):
+ """Return the ECSDataGranule elements as dictionary"""
+ self.getGranule()
+ datagran = {}
+ for i in self.granule.find('ECSDataGranule').getiterator():
+ if i.text.strip() != '':
+ datagran[i.tag] = i.text
+ return datagran
- def retRangeTime(self):
- """Return the RangeDateTime elements as dictionary
- """
- self.getGranule()
- rangeTime = {}
- for i in self.granule.find('RangeDateTime').getiterator():
- if i.text.strip() != '':
- rangeTime[i.tag] = i.text
- return rangeTime
+ def retPGEVersion(self):
+ """Return the PGEVersion element"""
+ self.getGranule()
+ return self.granule.find('PGEVersionClass').find('PGEVersion').text
- def retBoundary(self):
- """Return the maximum extend (Bounding Box) of the MODIS file as
- dictionary"""
- self.getGranule()
- self.boundary = []
- lat = []
- lon = []
- spatialContainer = self.granule.find('SpatialDomainContainer')
- horizontal = spatialContainer.find('HorizontalSpatialDomainContainer')
- boundary = horizontal.find('GPolygon').find('Boundary')
- for i in boundary.findall('Point'):
- la = float(i.find('PointLongitude').text)
- lo = float(i.find('PointLatitude').text)
- lon.append(la)
- lat.append(lo)
- self.boundary.append({'lat': la, 'lon': lo})
- extent = {'min_lat': min(lat), 'max_lat': max(lat), 'min_lon': min(lon),
- 'max_lon': max(lon)}
- return extent
+ def retRangeTime(self):
+ """Return the RangeDateTime elements as dictionary"""
+ self.getGranule()
+ rangeTime = {}
+ for i in self.granule.find('RangeDateTime').getiterator():
+ if i.text.strip() != '':
+ rangeTime[i.tag] = i.text
+ return rangeTime
- def retMeasure(self):
- """Return statistics of QA as dictionary"""
- value = {}
- self.getGranule()
- mes = self.granule.find('MeasuredParameter')
- mespc = mes.find('MeasuredParameterContainer')
- value['ParameterName'] = mespc.find('ParameterName').text
- meStat = mespc.find('QAStats')
- qastat = {}
- for i in meStat.getiterator():
- if i.tag != 'QAStats':
- qastat[i.tag] = i.text
- value['QAStats'] = qastat
- meFlag = mespc.find('QAFlags')
- flagstat = {}
- for i in meFlag.getiterator():
- if i.tag != 'QAFlags':
- flagstat[i.tag] = i.text
- value['QAFlags'] = flagstat
- return value
+ def retBoundary(self):
+ """Return the maximum extend (Bounding Box) of the MODIS file as
+ dictionary"""
+ self.getGranule()
+ self.boundary = []
+ lat = []
+ lon = []
+ spatialContainer = self.granule.find('SpatialDomainContainer')
+ horizontal = spatialContainer.find('HorizontalSpatialDomainContainer')
+ boundary = horizontal.find('GPolygon').find('Boundary')
+ for i in boundary.findall('Point'):
+ la = float(i.find('PointLongitude').text)
+ lo = float(i.find('PointLatitude').text)
+ lon.append(la)
+ lat.append(lo)
+ self.boundary.append({'lat': la, 'lon': lo})
+ extent = {'min_lat': min(lat), 'max_lat': max(lat),
+ 'min_lon': min(lon), 'max_lon': max(lon)}
+ return extent
- def retPlatform(self):
- """Return the platform values as dictionary."""
- value = {}
- self.getGranule()
- 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
- return value
+ def retMeasure(self):
+ """Return statistics of QA as dictionary"""
+ value = {}
+ self.getGranule()
+ mes = self.granule.find('MeasuredParameter')
+ mespc = mes.find('MeasuredParameterContainer')
+ value['ParameterName'] = mespc.find('ParameterName').text
+ meStat = mespc.find('QAStats')
+ qastat = {}
+ for i in meStat.getiterator():
+ if i.tag != 'QAStats':
+ qastat[i.tag] = i.text
+ value['QAStats'] = qastat
+ meFlag = mespc.find('QAFlags')
+ flagstat = {}
+ for i in meFlag.getiterator():
+ if i.tag != 'QAFlags':
+ flagstat[i.tag] = i.text
+ value['QAFlags'] = flagstat
+ return value
- def retPSA(self):
- """Return the PSA values as dictionary, the PSAName is the key and
- and PSAValue is the value
- """
- value = {}
- self.getGranule()
- psas = self.granule.find('PSAs')
- for i in psas.findall('PSA'):
- value[i.find('PSAName').text] = i.find('PSAValue').text
- return value
+ def retPlatform(self):
+ """Return the platform values as dictionary."""
+ value = {}
+ self.getGranule()
+ 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
+ return value
- def retInputGranule(self):
- """Return the input files (InputGranule) used to process the considered
- file"""
- value = []
- self.getGranule()
- for i in self.granule.find('InputGranule').getiterator():
- if i.tag != 'InputGranule':
- value.append(i.text)
- return value
+ def retPSA(self):
+ """Return the PSA values as dictionary, the PSAName is the key and
+ and PSAValue is the value
+ """
+ value = {}
+ self.getGranule()
+ psas = self.granule.find('PSAs')
+ for i in psas.findall('PSA'):
+ value[i.find('PSAName').text] = i.find('PSAValue').text
+ return value
- def retBrowseProduct(self):
- """Return the BrowseProduct element"""
- self.getGranule()
- try:
- value = self.granule.find('BrowseProduct').find('BrowseGranuleId').text
- except:
- value = None
- return value
+ def retInputGranule(self):
+ """Return the input files (InputGranule) used to process the considered
+ file"""
+ value = []
+ self.getGranule()
+ for i in self.granule.find('InputGranule').getiterator():
+ if i.tag != 'InputGranule':
+ value.append(i.text)
+ return value
- def confResample(self, spectral, res=None, output=None, datum='WGS84',
- resample='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 )',
- bound=None
- ):
- """Create the parameter file to use with resample MRT software to create
- tif (geotiff) file
+ def retBrowseProduct(self):
+ """Return the BrowseProduct element"""
+ self.getGranule()
+ try:
+ value = self.granule.find('BrowseProduct').find('BrowseGranuleId').text
+ except:
+ value = None
+ return value
- spectral = the spectral subset to be used, see the product table to
- understand the layer that you want use. For example:
+ def confResample(self, spectral, res=None, output=None, datum='WGS84',
+ resample='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 )', bound=None):
+ """Create the parameter file to use with resample MRT software to
+ create tif (geotiff) file
- - NDVI ( 1 1 1 0 0 0 0 0 0 0 0 0) copy only layer NDVI, EVI
- and QA VI the other layers are not used
- - LST ( 1 1 0 0 1 1 0 0 0 0 0 0 ) copy only layer daily and
- nightly temperature and QA
+ :param str spectral: the spectral subset to be used, see the product
+ table to understand the layer that you want use.
+ For example:
- res = the resolution for the output file, it must be set in the map
- unit of output projection system. The software will use the
- original resolution of input file if res not set
+ * NDVI ( 1 1 1 0 0 0 0 0 0 0 0 0) copy only layer
+ NDVI, EVI and QA VI the other layers are not used
+ * LST ( 1 1 0 0 1 1 0 0 0 0 0 0 ) copy only layer
+ daily and nightly temperature and QA
- output = the output name, if not set if not set the prefix name
- of input hdf file will be used
+ :param int res: the resolution for the output file, it must be set in
+ the map unit of output projection system. The software
+ will use the original resolution of input file if res
+ not set
- utm = the UTM zone if projection system is UTM
+ :param str output: the output name, if not set if not set the prefix
+ name of input hdf file will be used
- resample = the type of resampling, the valid values are:
- - NN (nearest neighbor)
- - BI (bilinear)
- - CC (cubic convolution)
+ :param utm: the UTM zone if projection system is UTM
- projtype = the output projection system, the valid values are:
- - AEA (Albers Equal Area)
- - ER (Equirectangular)
- - GEO (Geographic Latitude/Longitude)
- - HAM (Hammer)
- - ISIN (Integerized Sinusoidal)
- - IGH (Interrupted Goode Homolosine)
- - LA (Lambert Azimuthal)
- - LCC (LambertConformal Conic)
- - MERCAT (Mercator)
- - MOL (Mollweide)
- - PS (Polar Stereographic)
- - SIN (Sinusoidal)
- - UTM (Universal TransverseMercator)
+ :param str resample: the type of resampling, the valid values are:
- datum = the datum to use, the valid values are:
- - NAD27
- - NAD83
- - WGS66
- - WGS76
- - WGS84
- - NODATUM
+ * NN (nearest neighbor)
+ * BI (bilinear)
+ * CC (cubic convolution)
- projpar = a list of projection parameters, for more info check the
- Appendix C of MODIS reprojection tool user manual
- https://lpdaac.usgs.gov/content/download/4831/22895/file/mrt41_usermanual_032811.pdf
+ :param str projtype: the output projection system, valid values are:
- bound = dictionary with the following keys:
- - max_lat
- - max_lon
- - min_lat
- - min_lon
+ * AEA (Albers Equal Area)
+ * ER (Equirectangular)
+ * GEO (Geographic Latitude/Longitude)
+ * HAM (Hammer)
+ * ISIN (Integerized Sinusoidal)
+ * IGH (Interrupted Goode Homolosine)
+ * LA (Lambert Azimuthal)
+ * LCC (LambertConformal Conic)
+ * MERCAT (Mercator)
+ * MOL (Mollweide)
+ * PS (Polar Stereographic)
+ * SIN (Sinusoidal)
+ * UTM (Universal TransverseMercator)
+
+ :param str datum: the datum to use, the valid values are:
+
+ * NAD27
+ * NAD83
+ * WGS66
+ * WGS76
+ * WGS84
+ * NODATUM
+
+ :param str projpar: a list of projection parameters, for more info
+ check the Appendix C of MODIS reprojection tool
+ user manual https://lpdaac.usgs.gov/content/download/4831/22895/file/mrt41_usermanual_032811.pdf
+
+ :param dict bound: dictionary with the following keys:
+
+ * max_lat
+ * max_lon
+ * min_lat
+ * min_lon
"""
- #check if spectral it's write with correct construct ( value )
- if string.find(spectral, '(') == -1 or string.find(spectral, ')') == -1:
- raise IOError('ERROR: The spectral string should be similar to: ( 1 0 )')
- # output name
- if not output:
- fileout = self.tifname
- else:
- fileout = output
- # the name of the output parameters files for resample MRT software
- filename = os.path.join(self.path, '%s_mrt_resample.conf' % self.code)
- # if the file already exists it remove it
- if os.path.exists(filename):
- os.remove(filename)
- # open the file
- conFile = open(filename, 'w')
- conFile.write("INPUT_FILENAME = %s\n" % self.hdfname)
- conFile.write("SPECTRAL_SUBSET = %s\n" % spectral)
- conFile.write("SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG\n")
- if not bound:
- # return the boundary from the input xml file
- bound = self.retBoundary()
- else:
- if 'max_lat' not in bound or 'min_lat' not in bound or \
- 'min_lon' not in bound or 'max_lon' not in bound:
- raise IOError('bound variable is a dictionary with the following ' \
- 'keys: max_lat, min_lat, min_lon, max_lon')
- # Order: UL: N W - LR: S E
- conFile.write("SPATIAL_SUBSET_UL_CORNER = ( %f %f )\n" % (bound['max_lat'],
- bound['min_lon']))
- conFile.write("SPATIAL_SUBSET_LR_CORNER = ( %f %f )\n" % (bound['min_lat'],
- bound['max_lon']))
- conFile.write("OUTPUT_FILENAME = %s\n" % fileout)
- # if resample is in resam_list set the parameter otherwise return an error
- if resample in RESAM_LIST:
- conFile.write("RESAMPLING_TYPE = %s\n" % resample)
- else:
- raise IOError('The resampling type %s is not supportet.\n' \
- 'The resampling type supported are %s' % (resample,
- RESAM_LIST))
- # if projtype is in proj_list set the parameter otherwise return an error
- if projtype in PROJ_LIST:
- conFile.write("OUTPUT_PROJECTION_TYPE = %s\n" % projtype)
- else:
- raise IOError('The projection type %s is not supported.\n' \
- 'The projections supported are %s' % (projtype, PROJ_LIST))
- conFile.write("OUTPUT_PROJECTION_PARAMETERS = %s\n" % projpar)
- # if datum is in datum_list set the parameter otherwise return an error
- if datum in DATUM_LIST:
- conFile.write("DATUM = %s\n" % datum)
- else:
- raise IOError('The datum %s is not supported.\n' \
- 'The datum supported are %s' % (datum, DATUM_LIST))
- # if utm is not None write the UTM_ZONE parameter in the file
- if utm:
- conFile.write("UTM_ZONE = %s\n" % utm)
- # if res is not None write the OUTPUT_PIXEL_SIZE parameter in the file
- if res:
- conFile.write("OUTPUT_PIXEL_SIZE = %i\n" % res)
- conFile.close()
- return filename
+ # check if spectral it's write with correct construct ( value )
+ if string.find(spectral, '(') == -1 or string.find(spectral, ')') == -1:
+ raise Exception('ERROR: The spectral string should be similar to:'
+ ' ( 1 0 )')
+ # output name
+ if not output:
+ fileout = self.tifname
+ else:
+ fileout = output
+ # the name of the output parameters files for resample MRT software
+ filename = os.path.join(self.path, '{cod}_mrt_resample.conf'.format(cod=self.code))
+ # if the file already exists it remove it
+ if os.path.exists(filename):
+ os.remove(filename)
+ # open the file
+ conFile = open(filename, 'w')
+ conFile.write("INPUT_FILENAME = {name}\n".format(name=self.hdfname))
+ conFile.write("SPECTRAL_SUBSET = {spec}\n".format(spec=spectral))
+ conFile.write("SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG\n")
+ if not bound:
+ # return the boundary from the input xml file
+ bound = self.retBoundary()
+ else:
+ if 'max_lat' not in bound or 'min_lat' not in bound or \
+ 'min_lon' not in bound or 'max_lon' not in bound:
+ raise Exception('bound variable is a dictionary with the '
+ 'following keys: max_lat, min_lat, min_lon,'
+ ' max_lon')
+ # Order: UL: N W - LR: S E
+ conFile.write("SPATIAL_SUBSET_UL_CORNER = ( {mala} {milo} )"
+ "\n".format(mala=bound['max_lat'], milo=bound['min_lon']))
+ conFile.write("SPATIAL_SUBSET_LR_CORNER = ( {mila} {malo} )"
+ "\n".format(mila=bound['min_lat'], malo=bound['max_lon']))
+ conFile.write("OUTPUT_FILENAME = {out}\n".format(out=fileout))
+ # if resample is in resam_list set it otherwise return an error
+ if resample in RESAM_LIST:
+ conFile.write("RESAMPLING_TYPE = {res}\n".format(res=resample))
+ else:
+ raise Exception('The resampling type {res} is not supportet.\n'
+ 'The resampling type supported are '
+ '{reslist}'.format(res=resample, reslist=RESAM_LIST))
+ # if projtype is in proj_list set it otherwise return an error
+ if projtype in PROJ_LIST:
+ conFile.write("OUTPUT_PROJECTION_TYPE = {typ}\n".format(typ=projtype))
+ else:
+ raise Exception('The projection type {typ} is not supported.\n'
+ 'The projections supported are '
+ '{proj}'.format(typ=projtype, proj=PROJ_LIST))
+ conFile.write("OUTPUT_PROJECTION_PARAMETERS = {proj}\n".format(proj=projpar))
+ # if datum is in datum_list set the parameter otherwise return an error
+ if datum in DATUM_LIST:
+ conFile.write("DATUM = {dat}\n".format(dat=datum))
+ else:
+ raise Exception('The datum {dat} is not supported.\n'
+ 'The datum supported are '
+ '{datum}'.format(dat=datum, datum=DATUM_LIST))
+ # if utm is not None write the UTM_ZONE parameter in the file
+ if utm:
+ conFile.write("UTM_ZONE = {zone}\n".format(zone=utm))
+ # if res is not None write the OUTPUT_PIXEL_SIZE parameter in the file
+ if res:
+ conFile.write("OUTPUT_PIXEL_SIZE = {pix}\n".format(pix=res))
+ conFile.close()
+ return filename
- def confResample_swath(self, sds, geoloc, res, output=None,
- sphere='8', resample='NN', 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',
- bound=None
- ):
- """Create the parameter file to use with resample MRT software to create
- tif (geotiff) file
+ def confResample_swath(self, sds, geoloc, res, output=None,
+ sphere='8', resample='NN', 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', bound=None):
+ """Create the parameter file to use with resample MRT software to
+ create tif (geotiff) file
- sds = Name of band/s (Science Data Set) to resample
+ :param str sds: Name of band/s (Science Data Set) to resample
+ :param str geoloc: Name geolocation file (example MOD3, MYD3)
+ :param int res: the resolution for the output file, it must be set in
+ the map unit of output projection system. The software
+ will use the original resolution of input file if res
+ not set
- geoloc = Name geolocation file (example MOD3, MYD3)
+ :param str output: the output name, if not set the prefix name of
+ input hdf file will be used
- res = the resolution for the output file, it must be set in the map
- unit of output projection system. The software will use the
- original resolution of input file if res not set
+ :param int sphere: Output sphere number. Valid options are:
- output = the output name, if not set the prefix name
- of input hdf file will be used
+ * 0=Clarke 1866
+ * 1=Clarke 1880
+ * 2=Bessel
+ * 3=International 1967
+ * 4=International 1909
+ * 5=WGS 72
+ * 6=Everest
+ * 7=WGS 66
+ * 8=GRS1980/WGS 84
+ * 9=Airy
+ * 10=Modified Everest
+ * 11=Modified Airy
+ * 12=Walbeck
+ * 13=Southeast Asia
+ * 14=Australian National
+ * 15=Krassovsky
+ * 16=Hough
+ * 17=Mercury1960
+ * 18=Modified Mercury1968
+ * 19=Sphere 19 (Radius 6370997)
+ * 20=MODIS Sphere (Radius 6371007.181)
- sphere = Output sphere number. Valid options are:
- - 0=Clarke 1866
- - 1=Clarke 1880
- - 2=Bessel
- - 3=International 1967
- - 4=International 1909
- - 5=WGS 72
- - 6=Everest
- - 7=WGS 66
- - 8=GRS1980/WGS 84
- - 9=Airy
- - 10=Modified Everest
- - 11=Modified Airy
- - 12=Walbeck
- - 13=Southeast Asia
- - 14=Australian National
- - 15=Krassovsky
- - 16=Hough
- - 17=Mercury1960
- - 18=Modified Mercury1968
- - 19=Sphere 19 (Radius 6370997)
- - 20=MODIS Sphere (Radius 6371007.181)
+ :param str resample: the type of resampling, the valid values are:
- resample = the type of resampling, the valid values are:
- - NN (nearest neighbor)
- - BI (bilinear)
- - CC (cubic convolution)
+ * NN (nearest neighbor)
+ * BI (bilinear)
+ * CC (cubic convolution)
- projtype = the output projection system, the valid values are:
- - AEA (Albers Equal Area)
- - ER (Equirectangular)
- - GEO (Geographic Latitude/Longitude)
- - HAM (Hammer)
- - ISIN (Integerized Sinusoidal)
- - IGH (Interrupted Goode Homolosine)
- - LA (Lambert Azimuthal)
- - LCC (LambertConformal Conic)
- - MERCAT (Mercator)
- - MOL (Mollweide)
- - PS (Polar Stereographic),
- - SIN ()Sinusoidal)
- - UTM (Universal TransverseMercator)
+ :param str projtype: the output projection system, valid values are:
- utm = the UTM zone if projection system is UTM
+ * AEA (Albers Equal Area)
+ * ER (Equirectangular)
+ * GEO (Geographic Latitude/Longitude)
+ * HAM (Hammer)
+ * ISIN (Integerized Sinusoidal)
+ * IGH (Interrupted Goode Homolosine)
+ * LA (Lambert Azimuthal)
+ * LCC (LambertConformal Conic)
+ * MERCAT (Mercator)
+ * MOL (Mollweide)
+ * PS (Polar Stereographic),
+ * SIN ()Sinusoidal)
+ * UTM (Universal TransverseMercator)
- projpar = a list of projection parameters, for more info check
- the Appendix C of MODIS reprojection tool user manual
- https://lpdaac.usgs.gov/content/download/4831/22895/file/mrt41_usermanual_032811.pdf
+ :param utm: the UTM zone if projection system is UTM
- bound = dictionary with the following keys:
- - max_lat
- - max_lon
- - min_lat
- - min_lon
+ :param str projpar: a list of projection parameters, for more info
+ check the Appendix C of MODIS reprojection tool
+ user manual https://lpdaac.usgs.gov/content/download/4831/22895/file/mrt41_usermanual_032811.pdf
+
+ :param dict bound: dictionary with the following keys:
+
+ * max_lat
+ * max_lon
+ * min_lat
+ * min_lon
"""
- # output name
- if not output:
- fileout = self.tifname
- else:
- fileout = output
- # the name of the output parameters files for resample MRT software
- filename = os.path.join(self.path, '%s_mrt_resample.prm' % self.code)
- # if the file already exists it remove it
- if os.path.exists(filename):
- os.remove(filename)
- # open the file
- conFile = open(filename, 'w')
- conFile.write("INPUT_FILENAME = %s\n" % self.hdfname)
- conFile.write("GEOLOCATION_FILENAME = %s\n" % geoloc)
- conFile.write("INPUT_SDS_NAME = %s\n" % sds)
- conFile.write("OUTPUT_SPATIAL_SUBSET_TYPE = LAT_LONG\n")
- if not bound:
- # return the boundary from the input xml file
- bound = self.retBoundary()
- else:
- if 'max_lat' not in bound or 'min_lat' not in bound or \
- 'min_lon' not in bound or 'max_lon' not in bound:
- raise IOError('bound variable is a dictionary with the following ' \
- 'keys: max_lat, min_lat, min_lon, max_lon')
- # Order: UL: N W - LR: S E
- conFile.write("OUTPUT_SPACE_UPPER_LEFT_CORNER (LONG LAT) = %f %f\n" % (bound['max_lat'],
- bound['min_lon']))
- conFile.write("OUTPUT_SPACE_LOWER_RIGHT_CORNER (LONG LAT) = %f %f\n" % (bound['min_lat'],
- bound['max_lon']))
- conFile.write("OUTPUT_FILENAME = %s\n" % fileout)
- conFile.write("OUTPUT_FILE_FORMAT = GEOTIFF_FMT\n")
- # if resample is in resam_list set the parameter otherwise return an error
- if resample in RESAM_LIST_SWATH:
- conFile.write("KERNEL_TYPE (CC/BI/NN) = %s\n" % resample)
- else:
- raise IOError('The resampling type %s is not supportet.\n' \
- 'The resampling type supported are %s' % (resample,
- RESAM_LIST_SWATH))
- # if projtype is in proj_list set the parameter otherwise return an error
- if projtype in PROJ_LIST:
- conFile.write("OUTPUT_PROJECTION_NUMBER = %s\n" % projtype)
- else:
- raise IOError('The projection type %s is not supported.\n' \
- 'The projections supported are %s' % (projtype, PROJ_LIST))
- conFile.write("OUTPUT_PROJECTION_PARAMETER = %s\n" % projpar)
- # if sphere is in sphere_list set the parameter otherwise return an error
- if int(sphere) in SPHERE_LIST:
- conFile.write("OUTPUT_PROJECTION_SPHERE = %s\n" % sphere)
- else:
- raise IOError('The sphere %s is not supported.\n' \
- 'The spheres supported are %s' % (sphere, SPHERE_LIST))
- # if utm is not None write the UTM_ZONE parameter in the file
- if utm:
- if utm < '-60' or utm > '60':
- raise IOError('The valid UTM zone are -60 to 60')
- else:
- conFile.write("OUTPUT_PROJECTION_ZONE = %s\n" % utm)
- # if res is not None write the OUTPUT_PIXEL_SIZE parameter in the file
- if res:
- conFile.write("OUTPUT_PIXEL_SIZE = %f\n" % res)
- conFile.close()
- return filename
+ # output name
+ if not output:
+ fileout = self.tifname
+ else:
+ fileout = output
+ # the name of the output parameters files for resample MRT software
+ filename = os.path.join(self.path,
+ '{cod}_mrt_resample.prm'.format(cod=self.code))
+ # if the file already exists it remove it
+ if os.path.exists(filename):
+ os.remove(filename)
+ # open the file
+ conFile = open(filename, 'w')
+ conFile.write("INPUT_FILENAME = {name}\n".format(name=self.hdfname))
+ conFile.write("GEOLOCATION_FILENAME = {name}\n".format(name=geoloc))
+ conFile.write("INPUT_SDS_NAME = {name}\n".format(name=sds))
+ conFile.write("OUTPUT_SPATIAL_SUBSET_TYPE = LAT_LONG\n")
+ if not bound:
+ # return the boundary from the input xml file
+ bound = self.retBoundary()
+ else:
+ if 'max_lat' not in bound or 'min_lat' not in bound or \
+ 'min_lon' not in bound or 'max_lon' not in bound:
+ raise Exception('bound variable is a dictionary with the '
+ 'following keys: max_lat, min_lat, min_lon,'
+ ' max_lon')
+ # Order: UL: N W - LR: S E
+ conFile.write("OUTPUT_SPACE_UPPER_LEFT_CORNER (LONG LAT) = {milo} "
+ "{mala}\n".format(mala=bound['max_lat'],
+ milo=bound['min_lon']))
+ conFile.write("OUTPUT_SPACE_LOWER_RIGHT_CORNER (LONG LAT) = {mila} "
+ "{malo}\n".format(mila=bound['min_lat'],
+ malo=bound['max_lon']))
+ conFile.write("OUTPUT_FILENAME = {name}\n".format(name=fileout))
+ conFile.write("OUTPUT_FILE_FORMAT = GEOTIFF_FMT\n")
+ # if resample is in resam_list set it otherwise return an error
+ if resample in RESAM_LIST_SWATH:
+ conFile.write("KERNEL_TYPE (CC/BI/NN) = {res}"
+ "\n".format(res=resample))
+ else:
+ raise Exception('The resampling type {typ} is not supportet.\n'
+ 'The resampling type supported are '
+ '{swa}'.format(typ=resample, swa=RESAM_LIST_SWATH))
+ # if projtype is in proj_list set it otherwise return an error
+ if projtype in PROJ_LIST:
+ conFile.write("OUTPUT_PROJECTION_NUMBER = {typ}\n".format(typ=projtype))
+ else:
+ raise Exception('The projection type {typ} is not supported.\n'
+ 'The projections supported are '
+ '{proj}'.format(typ=projtype, proj=PROJ_LIST))
+ conFile.write("OUTPUT_PROJECTION_PARAMETER = {prj}\n".format(prj=projpar))
+ # if sphere is in sphere_list set it otherwise return an error
+ if int(sphere) in SPHERE_LIST:
+ conFile.write("OUTPUT_PROJECTION_SPHERE = {sph}\n".format(sph=sphere))
+ else:
+ raise Exception('The sphere {sph} is not supported.\nThe spheres'
+ 'supported are {sphere}'.format(sph=sphere,
+ sphere=SPHERE_LIST))
+ # if utm is not None write the UTM_ZONE parameter in the file
+ if utm:
+ if utm < '-60' or utm > '60':
+ raise Exception('The valid UTM zone are -60 to 60')
+ else:
+ conFile.write("OUTPUT_PROJECTION_ZONE = {utm}\n".format(utm=utm))
+ # if res is not None write the OUTPUT_PIXEL_SIZE parameter in the file
+ if res:
+ conFile.write("OUTPUT_PIXEL_SIZE = {res}\n".format(res=res))
+ conFile.close()
+ return filename
class parseModisMulti:
- """A class to obtain some variables for the xml file of several MODIS tiles.
- It can also create the xml file
- """
+ """A class to obtain some variables for the xml file of several MODIS
+ tiles. It can also create the xml file
- def __init__(self, hdflist):
- """hdflist = python list containing the hdf files"""
- from xml.etree import ElementTree
- self.ElementTree = ElementTree
- self.hdflist = hdflist
- self.parModis = []
- self.nfiles = 0
- # for each hdf files create a parseModis object
- for i in hdflist:
- self.parModis.append(parseModis(i))
- self.nfiles += 1
+ :param list hdflist: python list containing the hdf files
+ """
+ def __init__(self, hdflist):
+ """Function to initialize the object"""
+ from xml.etree import ElementTree
+ self.ElementTree = ElementTree
+ self.hdflist = hdflist
+ self.parModis = []
+ self.nfiles = 0
+ # for each hdf files create a parseModis object
+ for i in hdflist:
+ self.parModis.append(parseModis(i))
+ self.nfiles += 1
- def _most_common(self, lst):
- """Return the most common value of a list"""
- return max(set(lst), key=lst.count)
+ def _most_common(self, lst):
+ """Return the most common value of a list"""
+ return max(set(lst), key=lst.count)
+ def _checkval(self, vals):
+ """Internal function to return values from list
- def _checkval(self, vals):
- """Internal function to return values from list
+ :param list vals: list of values
+ """
+ if vals.count(vals[0]) == self.nfiles:
+ return [vals[0]]
+ else:
+ outvals = []
+ for i in vals:
+ if outvals.count(i) == 0:
+ outvals.append(i)
+ return outvals
- vals = list of values
- """
- if vals.count(vals[0]) == self.nfiles:
- return [vals[0]]
- else:
- outvals = []
- for i in vals:
- if outvals.count(i) == 0:
- outvals.append(i)
- return outvals
+ def _checkvaldict(self, vals):
+ """Internal function to return values from dictionary
- def _checkvaldict(self, vals):
- """Internal function to return values from dictionary
+ :param dict vals: dictionary of values
+ """
+ keys = vals[0].keys()
+ outvals = {}
+ for k in keys:
+ valtemp = []
+ for v in vals:
+ valtemp.append(v[k])
+ if valtemp.count(valtemp[0]) == self.nfiles:
+ outvals[k] = valtemp[0]
+ elif len(valtemp) == self.nfiles:
+ outvals[k] = self._most_common(valtemp)
+ else:
+ raise Exception('Something wrong reading XML files')
- vals = dictionary of values
- """
- keys = vals[0].keys()
- outvals = {}
- for k in keys:
- valtemp = []
- for v in vals:
- valtemp.append(v[k])
- if valtemp.count(valtemp[0]) == self.nfiles:
- outvals[k] = valtemp[0]
- elif len(valtemp) == self.nfiles:
- outvals[k] = self._most_common(valtemp)
- else:
- raise IOError('Something wrong reading XML files')
+ return outvals
- return outvals
+ def _minval(self, vals):
+ """Internal function to return the minimum value
- def _minval(self, vals):
- """Internal function to return the minimum value
+ :param list vals: list of values
+ """
+ outval = vals[0]
+ for i in range(1, len(vals)):
+ if outval > i:
+ outval = i
+ return outval
- vals = list of values
- """
- outval = vals[0]
- for i in range(1, len(vals)):
- if outval > i:
- outval = i
- return outval
+ def _maxval(self, vals):
+ """Internal function to return the maximum value
- def _maxval(self, vals):
- """Internal function to return the maximum value
+ :param list vals: list of values
+ """
+ outval = vals[0]
+ for i in range(1, len(vals)):
+ if outval < i:
+ outval = i
+ return outval
- vals = list of values
- """
- outval = vals[0]
- for i in range(1, len(vals)):
- if outval < i:
- outval = i
- return outval
+ def _cicle_values(self, obj, values):
+ """Internal function to add values from a dictionary
- def _cicle_values(self, obj, values):
- """Internal function to add values from a dictionary
+ :param obj: element to add values
- obj = element to add values
+ :param values: dictionary containing keys and values
+ """
+ for k, v in values.iteritems():
+ elem = self.ElementTree.SubElement(obj, k)
+ elem.text = v
- values = dictionary containing keys and values
- """
- for k,v in values.iteritems():
- elem = self.ElementTree.SubElement(obj, k)
- elem.text = v
+ def _addPoint(self, obj, lon, lat):
+ """Internal function to add a point in boundary xml tag
- def _addPoint(self, obj, lon, lat):
- """Internal function to add a point in boundary xml tag
+ :param obj: element to add point
+ :param lon: longitude of point
+ :param lat: latitude of point
+ """
+ pt = self.ElementTree.SubElement(obj, 'Point')
+ ptlon = self.ElementTree.SubElement(pt, 'PointLongitude')
+ ptlon.text = str(self.boundary[lon])
+ ptlat = self.ElementTree.SubElement(pt, 'PointLatitude')
+ ptlat.text = str(self.boundary[lat])
- obj = element to add point
+ def valDTD(self, obj):
+ """Function to add DTDVersion
- lon = longitude of point
+ :param obj: element to add DTDVersion
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retDTD())
+ for i in self._checkval(values):
+ dtd = self.ElementTree.SubElement(obj, 'DTDVersion')
+ dtd.text = i
- lat = latitude of point
- """
- pt = self.ElementTree.SubElement(obj, 'Point')
- ptlon = self.ElementTree.SubElement(pt, 'PointLongitude')
- ptlon.text = str(self.boundary[lon])
- ptlat = self.ElementTree.SubElement(pt, 'PointLatitude')
- ptlat.text = str(self.boundary[lat])
+ def valDataCenter(self, obj):
+ """Function to add DataCenter
- def valDTD(self, obj):
- """Function to add DTDVersion
+ :param obj: element to add DataCenter
+ """
+ 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
- obj = element to add DTDVersion
- """
- values = []
- for i in self.parModis:
- values.append(i.retDTD())
- for i in self._checkval(values):
- dtd = self.ElementTree.SubElement(obj, 'DTDVersion')
- dtd.text = i
+ def valGranuleUR(self, obj):
+ """Function to add GranuleUR
- def valDataCenter(self, obj):
- """Function to add DataCenter
+ :param obj: element to add GranuleUR
+ """
+ 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
- obj = element to add DataCenter
- """
- 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 valDbID(self, obj):
+ """Function to add DbID
- def valGranuleUR(self, obj):
- """Function to add GranuleUR
+ :param obj: element to add DbID
+ """
+ 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
- obj = element to add GranuleUR
- """
- 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 valInsTime(self, obj):
+ """Function to add the minimum of InsertTime
- def valDbID(self, obj):
- """Function to add DbID
+ :param obj: element to add InsertTime
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retInsertTime())
+ obj.text = self._minval(values)
- obj = element to add DbID
- """
- 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 valCollectionMetaData(self, obj):
+ """Function to add CollectionMetaData
- def valInsTime(self, obj):
- """Function to add the minimum of InsertTime
+ :param obj: element to add CollectionMetaData
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retCollectionMetaData())
+ self._cicle_values(obj, self._checkvaldict(values))
- obj = element to add InsertTime
- """
- values = []
- for i in self.parModis:
- values.append(i.retInsertTime())
- obj.text = self._minval(values)
+ def valDataFiles(self, obj):
+ """Function to add DataFileContainer
- def valCollectionMetaData(self, obj):
- """Function to add CollectionMetaData
+ :param obj: element to add DataFileContainer
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retDataFiles())
+ for i in values:
+ dfc = self.ElementTree.SubElement(obj, 'DataFileContainer')
+ self._cicle_values(dfc, i)
- obj = element to add CollectionMetaData
- """
- values = []
- for i in self.parModis:
- values.append(i.retCollectionMetaData())
- self._cicle_values(obj, self._checkvaldict(values))
+ def valPGEVersion(self, obj):
+ """Function to add PGEVersion
- def valDataFiles(self, obj):
- """Function to add DataFileContainer
+ :param obj: element to add PGEVersion
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retPGEVersion())
+ for i in self._checkval(values):
+ pge = self.ElementTree.SubElement(obj, 'PGEVersion')
+ pge.text = i
- obj = element to add DataFileContainer
- """
- values = []
- for i in self.parModis:
- values.append(i.retDataFiles())
- for i in values:
- dfc = self.ElementTree.SubElement(obj, 'DataFileContainer')
- self._cicle_values(dfc, i)
+ def valRangeTime(self, obj):
+ """Function to add RangeDateTime
- def valPGEVersion(self, obj):
- """Function to add PGEVersion
-
- obj = element to add PGEVersion
- """
- values = []
- for i in self.parModis:
- values.append(i.retPGEVersion())
- for i in self._checkval(values):
- pge = self.ElementTree.SubElement(obj, 'PGEVersion')
- pge.text = i
+ :param obj: element to add RangeDateTime
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retRangeTime())
+ self._cicle_values(obj, self._checkvaldict(values))
- def valRangeTime(self, obj):
- """Function to add RangeDateTime
+ def valBound(self):
+ """Function return the Bounding Box of mosaic"""
+ boundary = self.parModis[0].retBoundary()
+ for i in range(1, len(self.parModis)):
+ bound = self.parModis[i].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']
+ self.boundary = boundary
- obj = element to add RangeDateTime
- """
- values = []
- for i in self.parModis:
- values.append(i.retRangeTime())
- self._cicle_values(obj, self._checkvaldict(values))
+ def valMeasuredParameter(self, obj):
+ """Function to add ParameterName
- def valBound(self):
- """Function return the Bounding Box of mosaic
- """
- boundary = self.parModis[0].retBoundary()
- for i in range(1, len(self.parModis)):
- bound = self.parModis[i].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']
- self.boundary = boundary
+ :param obj: element to add ParameterName
+ """
+ valuesQAStats = []
+ valuesQAFlags = []
+ valuesParameter = []
+ for i in self.parModis:
+ valuesQAStats.append(i.retMeasure()['QAStats'])
+ valuesQAFlags.append(i.retMeasure()['QAFlags'])
+ valuesParameter.append(i.retMeasure()['ParameterName'])
+ for i in self._checkval(valuesParameter):
+ pn = self.ElementTree.SubElement(obj, 'ParameterName')
+ pn.text = i
- def valMeasuredParameter(self, obj):
- """Function to add ParameterName
+ def valInputPointer(self, obj):
+ """Function to add InputPointer
- obj = element to add ParameterName
- """
- valuesQAStats = []
- valuesQAFlags = []
- valuesParameter = []
- for i in self.parModis:
- valuesQAStats.append(i.retMeasure()['QAStats'])
- valuesQAFlags.append(i.retMeasure()['QAFlags'])
- valuesParameter.append(i.retMeasure()['ParameterName'])
- for i in self._checkval(valuesParameter):
- pn = self.ElementTree.SubElement(obj, 'ParameterName')
- pn.text = i
+ :param obj: element to add InputPointer
+ """
+ for i in self.parModis:
+ for v in i.retInputGranule():
+ ip = self.ElementTree.SubElement(obj, 'InputPointer')
+ ip.text = v
- def valInputPointer(self, obj):
- """Function to add InputPointer
+ def valPlatform(self, obj):
+ """Function to add Platform elements
- obj = element to add InputPointer
- """
- for i in self.parModis:
- for v in i.retInputGranule():
- ip = self.ElementTree.SubElement(obj, 'InputPointer')
- ip.text = v
+ :param obj: element to add Platform elements
+ """
+ valuesSName = []
+ valuesInstr = []
+ valuesSensor = []
+ for i in self.parModis:
+ valuesSName.append(i.retPlatform()['PlatformShortName'])
+ valuesInstr.append(i.retPlatform()['InstrumentShortName'])
+ valuesSensor.append(i.retPlatform()['SensorShortName'])
+ for i in self._checkval(valuesSName):
+ pn = self.ElementTree.SubElement(obj, 'PlatformShortName')
+ pn.text = i
- def valPlatform(self, obj):
- """Function to add Platform elements
+ valInstr = self._checkval(valuesInstr)
+ valSens = self._checkval(valuesSensor)
- obj = element to add Platform elements
- """
- valuesSName = []
- valuesInstr = []
- valuesSensor = []
- for i in self.parModis:
- valuesSName.append(i.retPlatform()['PlatformShortName'])
- valuesInstr.append(i.retPlatform()['InstrumentShortName'])
- valuesSensor.append(i.retPlatform()['SensorShortName'])
- for i in self._checkval(valuesSName):
- pn = self.ElementTree.SubElement(obj, 'PlatformShortName')
- pn.text = i
+ if len(valInstr) != len(valSens):
+ raise Exception('Something wrong reading XML files')
+ else:
+ for i in range(len(valInstr)):
+ ins = self.ElementTree.SubElement(obj, 'Instrument')
+ pn = self.ElementTree.SubElement(ins, 'InstrumentShortName')
+ pn.text = valInstr[i]
+ sens = self.ElementTree.SubElement(ins, 'Sensor')
+ ps = self.ElementTree.SubElement(sens, 'SensorShortName')
+ ps.text = valSens[i]
- valInstr = self._checkval(valuesInstr)
- valSens = self._checkval(valuesSensor)
+ def valInsertTime(self, obj):
+ """Function to add InsertTime elements
- if len(valInstr) != len(valSens):
- raise IOError('Something wrong reading XML files')
- else:
- for i in range(len(valInstr)):
- ins = self.ElementTree.SubElement(obj, 'Instrument')
- pn = self.ElementTree.SubElement(ins, 'InstrumentShortName')
- pn.text = valInstr[i]
- sens = self.ElementTree.SubElement(ins, 'Sensor')
- ps = self.ElementTree.SubElement(sens, 'SensorShortName')
- ps.text = valSens[i]
+ :param obj: element to add InsertTime elements
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retInsertTime())
+ for i in self._checkval(values):
+ gur = self.ElementTree.SubElement(obj, 'InsertTime')
+ gur.text = i
- def writexml(self, outputname):
- """Write a xml file for a mosaic
+ def valLastUpdate(self, obj):
+ """Function to add LastUpdate elements
- outputname = the name of xml file
- """
- # the root element
- granule = self.ElementTree.Element('GranuleMetaDataFile')
- # add DTDVersion
- self.valDTD(granule)
- # add DataCenterId
- self.valDataCenter(granule)
- # add GranuleURMetaData
- gurmd = self.ElementTree.SubElement(granule, 'GranuleURMetaData')
- # add GranuleUR
- self.valGranuleUR(gurmd)
- # add dbID
- self.valDbID(gurmd)
+ :param obj: element to add LastUpdate elements
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retLastUpdate())
+ for i in self._checkval(values):
+ gur = self.ElementTree.SubElement(obj, 'LastUpdate')
+ gur.text = i
- # TODO ADD InsertTime LastUpdate
+ def valDataGranule(self, obj):
+ """Function to add DataFileContainer
- # add CollectionMetaData
- cmd = self.ElementTree.SubElement(gurmd, 'CollectionMetaData')
- self.valCollectionMetaData(cmd)
- # add DataFiles
- df = self.ElementTree.SubElement(gurmd, 'DataFiles')
- self.valDataFiles(df)
+ :param obj: element to add DataFileContainer
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retDataGranule())
+ for i in values:
+ dfc = self.ElementTree.SubElement(obj, 'ECSDataGranule')
+ self._cicle_values(dfc, i)
- # TODO ADD ECSDataGranule
+ def valBrowseProduct(self, obj):
+ """Function to add BrowseGranuleId
- # add PGEVersionClass
- pgevc = self.ElementTree.SubElement(gurmd, 'PGEVersionClass')
- self.valPGEVersion(pgevc)
- # add RangeDateTime
- rdt = self.ElementTree.SubElement(gurmd, 'RangeDateTime')
- self.valRangeTime(rdt)
- # SpatialDomainContainer
- sdc = self.ElementTree.SubElement(gurmd, 'SpatialDomainContainer')
- hsdc = self.ElementTree.SubElement(sdc, 'HorizontalSpatialDomainContainer')
- gp = self.ElementTree.SubElement(hsdc, 'GPolygon')
- bound = self.ElementTree.SubElement(gp, 'Boundary')
- self.valBound()
- self._addPoint(bound, 'min_lon', 'max_lat')
- self._addPoint(bound, 'max_lon', 'max_lat')
- self._addPoint(bound, 'min_lon', 'min_lat')
- self._addPoint(bound, 'max_lon', 'min_lat')
- # add MeasuredParameter
- mp = self.ElementTree.SubElement(gurmd, 'MeasuredParameter')
- mpc = self.ElementTree.SubElement(mp, 'MeasuredParameterContainer')
- self.valMeasuredParameter(mpc)
- # Platform
- pl = self.ElementTree.SubElement(gurmd, 'Platform')
- self.valPlatform(pl)
+ :param obj: element to add BrowseGranuleId
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retBrowseProduct())
+ for i in self._checkval(values):
+ dfc = self.ElementTree.SubElement(obj, 'BrowseGranuleId')
+ dfc.text = i
- # add PSAs
- psas = self.ElementTree.SubElement(gurmd, 'PSAs')
- # TODO ADD all PSA
+ def valPSA(self, obj):
+ """Function to add PSA
- # add InputGranule and InputPointer
- ig = self.ElementTree.SubElement(gurmd, 'InputGranule')
- self.valInputPointer(ig)
- # TODO ADD BrowseProduct
- output = open(outputname, '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(self.ElementTree.tostring(granule))
- output.close()
+ :param obj: element to add PSA
+ """
+ values = []
+ for i in self.parModis:
+ values.append(i.retPSA())
+ for k in sorted(values[0].keys()):
+ psa = self.ElementTree.SubElement(obj, 'PSA')
+ psaname = self.ElementTree.SubElement(psa, 'PSAName')
+ psaname.text = k
+ for s in values:
+ psaval = self.ElementTree.SubElement(psa, 'PSAValue')
+ psaval.text = s[k]
+
+ def writexml(self, outputname, pretty=True):
+ """Write a xml file for a mosaic
+
+ :param str outputname: the name of output xml file
+ :param bool pretty: write prettyfy output, by default true
+ """
+
+ # the root element
+ granule = self.ElementTree.Element('GranuleMetaDataFile')
+ # add DTDVersion
+ self.valDTD(granule)
+ # add DataCenterId
+ self.valDataCenter(granule)
+ # add GranuleURMetaData
+ gurmd = self.ElementTree.SubElement(granule, 'GranuleURMetaData')
+ # add GranuleUR
+ self.valGranuleUR(gurmd)
+ # add dbID
+ self.valDbID(gurmd)
+ # add InsertTime
+ self.valInsertTime(gurmd)
+ # add LastUpdate
+ self.valLastUpdate(gurmd)
+ # add CollectionMetaData
+ cmd = self.ElementTree.SubElement(gurmd, 'CollectionMetaData')
+ self.valCollectionMetaData(cmd)
+ # add DataFiles
+ df = self.ElementTree.SubElement(gurmd, 'DataFiles')
+ self.valDataFiles(df)
+ # add ECSDataGranule
+ self.valDataGranule(gurmd)
+ # add PGEVersionClass
+ pgevc = self.ElementTree.SubElement(gurmd, 'PGEVersionClass')
+ self.valPGEVersion(pgevc)
+ # add RangeDateTime
+ rdt = self.ElementTree.SubElement(gurmd, 'RangeDateTime')
+ self.valRangeTime(rdt)
+ # add SpatialDomainContainer
+ sdc = self.ElementTree.SubElement(gurmd, 'SpatialDomainContainer')
+ hsdc = self.ElementTree.SubElement(sdc, 'HorizontalSpatialDomainContainer')
+ gp = self.ElementTree.SubElement(hsdc, 'GPolygon')
+ bound = self.ElementTree.SubElement(gp, 'Boundary')
+ self.valBound()
+ self._addPoint(bound, 'min_lon', 'max_lat')
+ self._addPoint(bound, 'max_lon', 'max_lat')
+ self._addPoint(bound, 'min_lon', 'min_lat')
+ self._addPoint(bound, 'max_lon', 'min_lat')
+ # add MeasuredParameter
+ mp = self.ElementTree.SubElement(gurmd, 'MeasuredParameter')
+ mpc = self.ElementTree.SubElement(mp, 'MeasuredParameterContainer')
+ self.valMeasuredParameter(mpc)
+ # add Platform
+ pl = self.ElementTree.SubElement(gurmd, 'Platform')
+ self.valPlatform(pl)
+ # add PSAs
+ psas = self.ElementTree.SubElement(gurmd, 'PSAs')
+ # add all PSA
+ self.valPSA(psas)
+ # add InputGranule and InputPointer
+ ig = self.ElementTree.SubElement(gurmd, 'InputGranule')
+ self.valInputPointer(ig)
+ # add BrowseProduct
+ bp = self.ElementTree.SubElement(gurmd, 'BrowseProduct')
+ self.valBrowseProduct(bp)
+ output = open(outputname, '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">')
+ if pretty:
+ import xml.dom.minidom as minidom
+ reparsed = minidom.parseString(self.ElementTree.tostring(granule))
+ output.write(reparsed.toprettyxml(indent="\t"))
+ else:
+ output.write(self.ElementTree.tostring(granule))
+ output.close()
Modified: grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py
===================================================================
--- grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py 2015-10-19 01:03:34 UTC (rev 66525)
+++ grass-addons/grass7/raster/r.modis/libmodis/rmodislib.py 2015-10-19 13:43:18 UTC (rev 66526)
@@ -19,7 +19,7 @@
# interface to g.proj -p
-def get_proj():
+def get_proj(flag='p'):
"""!Returns the output from running "g.proj -p" plus towgs84 parameter
(g.proj -d), as a dictionary. Example:
@@ -31,17 +31,20 @@
@return dictionary of projection values
"""
- gproj = grass.read_command('g.proj', flags='p')
- listproj = gproj.split('\n')
- listproj.remove('-PROJ_INFO-------------------------------------------------')
- listproj.remove('-PROJ_UNITS------------------------------------------------')
- listproj.remove('')
- proj = {}
- for i in listproj:
- ilist = i.split(':')
- proj[ilist[0].strip()] = ilist[1].strip()
- proj.update(grass.parse_command('g.proj', flags='j'))
- return proj
+ gproj = grass.read_command('g.proj', flags=flag)
+ if flag == 'p':
+ listproj = gproj.split('\n')
+ listproj.remove('-PROJ_INFO-------------------------------------------------')
+ listproj.remove('-PROJ_UNITS------------------------------------------------')
+ listproj.remove('')
+ proj = {}
+ for i in listproj:
+ ilist = i.split(':')
+ proj[ilist[0].strip()] = ilist[1].strip()
+ proj.update(grass.parse_command('g.proj', flags='j'))
+ return proj
+ elif flag == 'w':
+ return gproj.replace('\n', '').replace(' ', '')
class product:
@@ -83,7 +86,7 @@
surf_specqa = '( 1 1 1 1 1 1 1 1 0 0 0 0 0 )'
surf_suff = {'.sur_refl_b01': '.sur_refl_qc_500m', '.sur_refl_b02':
'.sur_refl_qc_500m', '.sur_refl_b03': '.sur_refl_qc_500m',
- '.sur_refl_b04': '.sur_refl_qc_500m', '.sur_refl_b05':
+ '.sur_refl_b04': '.sur_refl_qc_500m', '.sur_refl_b05':
'.sur_refl_qc_500m', '.sur_refl_b06': '.sur_refl_qc_500m',
'.sur_refl_b07': '.sur_refl_qc_500m'}
Modified: grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.py
===================================================================
--- grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.py 2015-10-19 01:03:34 UTC (rev 66525)
+++ grass-addons/grass7/raster/r.modis/r.modis.download/r.modis.download.py 2015-10-19 13:43:18 UTC (rev 66526)
@@ -29,6 +29,10 @@
#% key: g
#% description: Return the name of file containing the list of HDF tiles downloaded in shell script style
#%end
+#%flag
+#% key: c
+#% description: Does not perform GDAL check on downloaded images
+#%end
#%option G_OPT_F_INPUT
#% key: settings
#% label: Full path to settings file
@@ -73,6 +77,7 @@
#% required: no
#%end
+
# import library
import os
import sys
@@ -116,10 +121,11 @@
int(firstSplit[2]))
lastSplit = second.split('-')
lastDay = date(int(lastSplit[0]), int(lastSplit[1]), int(lastSplit[2]))
+ if firstDay < lastDay:
+ grass.fatal(_("End day has to be bigger then start day"))
delta = firstDay - lastDay
valueDelta = int(delta.days)
return valueDay, second, valueDelta
-
# no set start and end day
if options['startday'] == '' and options['endday'] == '':
return None, None, 10
@@ -204,13 +210,17 @@
debug_opt = True
else:
debug_opt = False
+ if flags['c']:
+ checkgdal = False
+ else:
+ checkgdal = True
for produ in products:
prod = product(produ).returned()
#start modis class
-
modisOgg = downModis(url=prod['url'], user=user, password=passwd,
destinationFolder=fold, tiles=tiles, path=prod['folder'],
- today=firstday, enddate=finalday, delta=delta, debug=debug_opt)
+ today=firstday, enddate=finalday, delta=delta, debug=debug_opt,
+ checkgdal=checkgdal)
# connect to ftp
modisOgg.connect()
if modisOgg.nconnection <= 20:
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 2015-10-19 01:03:34 UTC (rev 66525)
+++ grass-addons/grass7/raster/r.modis/r.modis.import/r.modis.import.py 2015-10-19 13:43:18 UTC (rev 66526)
@@ -44,7 +44,7 @@
#% key_desc: path
#% description: Full path to MRT directory
#% gisprompt: old,dir,input
-#% required: yes
+#% required: no
#%end
#%option
#% key: dns
@@ -98,8 +98,9 @@
sys.path.append(path)
# try to import pymodis (modis) and some classes for r.modis.download
-from rmodislib import resampling, product, projection
-from convertmodis import convertModis, createMosaic
+from rmodislib import resampling, product, projection, get_proj
+from convertmodis import convertModis, createMosaic
+from convertmodis_gdal import createMosaicGDAL, convertModisGDAL
from parsemodis import parseModis
@@ -120,10 +121,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 is 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 is True:
day = line.split('/')[-1].split('.')[1]
if day in filelist:
filelist[day].append(line)
@@ -188,7 +189,7 @@
res = None
try:
conf = pm.confResample(spectr, res, pref, dat, resampl,
- proj, zone, projpar)
+ proj, zone, projpar)
return conf
except IOError, e:
grass.fatal(e)
@@ -256,8 +257,8 @@
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))
+ 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,
@@ -297,11 +298,26 @@
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)
- # create convertModis class and convert it in tif file
- execmodis = convertModis(hdf, confname, options['mrtpath'])
+ if options['mrtpath']:
+ # create conf file fro mrt tools
+ pm = parseModis(hdf)
+ confname = confile(pm, options, an)
+ # create convertModis class and convert it in tif file
+ execmodis = convertModis(hdf, confname, options['mrtpath'])
+ else:
+ projwkt = get_proj('w')
+ projObj = projection()
+ pref = listfile[0].split('/')[-1]
+ prod = product().fromcode(pref.split('.')[0])
+ spectr = spectral(options, prod, an)
+ if projObj.returned() != 'GEO':
+ res = int(prod['res']) * int(projObj.proj['meters'])
+ else:
+ res = None
+ prod = product().fromcode(pref.split('.')[0])
+ outname = "%s.%s.mosaic" % (pref.split('.')[0], pref.split('.')[1])
+ execmodis = convertModisGDAL(hdf, outname, spectr, res,
+ wkt=projwkt)
execmodis.run()
output = prefix(options)
if not output:
@@ -318,31 +334,47 @@
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.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
- prod = product().fromcode(listfiles[0].split('/')[-1].split('.')[0])
+ pref = listfiles[0].split('/')[-1]
+ prod = product().fromcode(pref.split('.')[0])
spectr = spectral(options, prod, an)
spectr = spectr.lstrip('( ').rstrip(' )')
+ outname = "%s.%s.mosaic" % (pref.split('.')[0], pref.split('.')[1])
# create mosaic
- cm = createMosaic(tempfile.name, outname, options['mrtpath'], spectr)
- cm.run()
+ if options['mrtpath']:
+ # create the file with the list of name
+ tempfile = open(os.path.join(targetdir, pid), 'w')
+ tempfile.writelines(listfiles)
+ tempfile.close()
+ # basedir of tempfile, where hdf files are write
+ basedir = os.path.split(tempfile.name)[0]
+ # return the spectral subset in according mrtmosaic tool format
+ cm = createMosaic(tempfile.name, outname, options['mrtpath'],
+ spectr)
+ cm.run()
+ else:
+ basedir = targetdir
+ cm = createMosaicGDAL(listfiles, spectr)
+ cm.write_vrt(outname)
# list of hdf files
hdfiles = glob.glob1(basedir, outname + "*.hdf")
for i in hdfiles:
# the full path to hdf file
hdf = os.path.join(basedir, i)
- # create conf file fro mrt tools
- pm = parseModis(hdf)
- confname = confile(pm, options, an, True)
# create convertModis class and convert it in tif file
- execmodis = convertModis(hdf, confname, options['mrtpath'])
+ if options['mrtpath']:
+ # create conf file fro mrt tools
+ pm = parseModis(hdf)
+ confname = confile(pm, options, an, True)
+ execmodis = convertModis(hdf, confname, options['mrtpath'])
+ else:
+ projwkt = get_proj('w')
+ projObj = projection()
+ if projwkt.returned() != 'GEO':
+ res = int(prod['res']) * int(projObj.proj['meters'])
+ else:
+ res = None
+ execmodis = convertModisGDAL(hdf, outname, res, wkt=projwkt,
+ vrt=True)
execmodis.run()
# remove hdf
if remove:
@@ -350,7 +382,7 @@
import_tif(outname, basedir, remove, ow, pm)
os.remove(hdf)
os.remove(hdf + '.xml')
- # or move the hdf and hdf.xml to the dir where are the original files
+ # move the hdf and hdf.xml to the dir where are the original files
else:
# import tif files
import_tif(outname, basedir, remove, ow, pm, targetdir)
@@ -373,14 +405,15 @@
return 0
# return an error if q and spectral are set
if not flags['q'] and options['spectral'] != '':
- grass.warning(_('If no QA layer chosen in the "spectral" option'\
- + ' the command will report an error'))
+ grass.warning(_('If no QA layer chosen in the "spectral" option'
+ ' the command will report an error'))
# return an error if both dns and files option are set or not
if options['dns'] == '' and options['files'] == '':
grass.fatal(_('Choose one of "dns" or "files" options'))
return 0
elif options['dns'] != '' and options['files'] != '':
- grass.fatal(_('It is not possible set "dns" and "files" options together'))
+ grass.fatal(_('It is not possible set "dns" and "files"'
+ ' options together'))
return 0
# check the version
version = grass.core.version()
@@ -404,7 +437,8 @@
analyze = 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'))
+ 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)
More information about the grass-commit
mailing list