Source code for pymodis.convertmodis_gdal

#!/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@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`

"""

# python 2 and 3 compatibility
from __future__ import print_function
from __future__ import division
from collections import OrderedDict
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]]'


[docs]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
[docs]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 preojection of output file :param str wkt: the WKT string for the preojection 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 isinstance(subset, list): self.subset = subset elif isinstance(subset, str): 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, src.GetProjection(), 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') 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') 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, quiet=False): """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 list(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 possible to create dataset %s' % out_name) 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, l_src_ds.GetProjection(), self.dst_wkt, self.resampling, 0, self.error_threshold, cbk, cbk_user_data) if not quiet: print("Layer {name} reprojected".format(name=l)) except: raise Exception('Not possible to reproject dataset ' '{name}'.format(name=l)) dst_ds.SetMetadata(meta) dst_ds = None l_src_ds = None return 0
[docs] def run_vrt_separated(self): """Reproject VRT created by createMosaicGDAL, function write_vrt with separated=True """ self._createWarped(self.in_name) self._reprojectOne(self.in_name) print("Dataset '{name}' reprojected".format(name=self.in_name))
[docs] def run(self, quiet=False): """Reproject all the subset of chosen 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], quiet=quiet) n = n + 1 if not quiet: print("All layer for dataset '{name}' " "reprojected".format(name=self.in_name))
# =============================================================================
[docs]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
[docs]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
[docs]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. """
[docs] 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 list(meta.keys()): self.fill_value = meta['_FillValue'] elif fh.GetRasterBand(1).GetNoDataValue(): self.fill_value = fh.GetRasterBand(1).GetNoDataValue() 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
[docs] 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)
[docs]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 isinstance(subset, list): self.subset = subset elif isinstance(subset, str): 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 chosen subset""" if isinstance(self.in_names, list): 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 = OrderedDict() 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] = list() 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 = OrderedDict() for k, v in self.layers.items(): 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 = list(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[list(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
[docs] 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)
[docs] def run(self, output, quiet=False): """Create the mosaic :param str output: the name of output file """ values = list(self.file_infos.values()) l1 = values[0][0] xsize, ysize, geotransform = self._calculateNewSize() t_fh = self.driver.Create(output, xsize, ysize, len(list(self.file_infos.keys())), l1.band_type) if t_fh is None: raise Exception('Not possible to create dataset %s' % output) t_fh.SetGeoTransform(geotransform) t_fh.SetProjection(l1.projection) i = 1 for names in list(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 if not quiet: print("The mosaic file {name} has been " "created".format(name=output)) return True
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
[docs] def write_vrt(self, output, separate=True, quiet=False): """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(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=gdal.GetDataTypeName(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 list(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(f, geot) out.write('\t</VRTRasterBand>\n') out.write('</VRTDataset>\n') out.close() else: values = list(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 list(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(f, geot) out.write('\t</VRTRasterBand>\n') band += 1 out.write('</VRTDataset>\n') out.close() if not quiet: print("The VRT mosaic file {name} has been " "created".format(name=output)) return True