# -*- coding: utf-8 -*-
"""**Functionality related to convert format file.**

InaSAFE Disaster risk assessment tool developed by AusAid and World Bank

Contact : ole.moller.nielsen@gmail.com

.. note:: 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.

Initially this was adapted from shake_event.py and now realtime uses this.
"""

__author__ = 'ismail@kartoza.com'
__version__ = '0.5.0'
__date__ = '11/02/2013'
__copyright__ = ('Copyright 2012, Australia Indonesia Facility for '
                 'Disaster Reduction')

import os
import sys
import shutil
import logging
import codecs
from xml.dom import minidom
from datetime import datetime
from pytz import timezone
from subprocess import call, CalledProcessError

from osgeo import gdal, ogr
from osgeo.gdalconst import GA_ReadOnly
# This import is required to enable PyQt API v2
# noinspection PyUnresolvedReferences
# pylint: disable=unused-import
import qgis
# pylint: enable=unused-import
from qgis.core import (
    QgsVectorLayer,
    QgsFeatureRequest,
    QgsRectangle)

from safe.common.utilities import which, romanise
from safe.common.exceptions import (
    GridXmlFileNotFoundError,
    GridXmlParseError,
    ContourCreationError,
    InvalidLayerError)
from safe.utilities.styling import mmi_colour
from safe.utilities.unicode import get_string

LOGGER = logging.getLogger('InaSAFE')


def data_dir():
    """Return the path to the standard data dir for e.g. geonames data

    :returns: Returns the default data directory.
    :rtype: str
    """
    dir_path = os.path.abspath(
        os.path.join(os.path.dirname(__file__), 'converter_data'))
    return dir_path


class ShakeGrid(object):
    """A converter for USGS shakemap grid.xml files to geotiff."""

    def __init__(
            self,
            title,
            source,
            grid_xml_path,
            output_dir=None,
            output_basename=None,
            algorithm_filename_flag=True):
        """Constructor.

        :param title: The title of the earthquake that will be also added to
            keywords file on some generated products.
        :type title: str

        :param source: The source of the earthquake that will be also added
            to keywords file on some generated products.
        :type source: str

        :param grid_xml_path: Path to grid XML file.
        :type grid_xml_path: str

        :param output_dir: mmi output directory
        :type output_dir: str

        :param output_basename: mmi file name without extension.
        :type output_basename: str

        :param algorithm_filename_flag: Flag whether to use the algorithm in
            the output file's name.
        :type algorithm_filename_flag: bool

        :returns: The instance of the class.
        :rtype: ShakeGrid

        :raises: EventXmlParseError
        """
        self.title = title
        self.source = source
        self.latitude = None
        self.longitude = None
        self.magnitude = None
        self.depth = None
        self.description = None
        self.location = None
        self.day = None
        self.month = None
        self.year = None
        self.time = None
        self.hour = None
        self.minute = None
        self.second = None
        self.time_zone = None
        self.x_minimum = None
        self.x_maximum = None
        self.y_minimum = None
        self.y_maximum = None
        # The bounding box of the grid as QgsRectangle
        self.grid_bounding_box = None
        self.rows = None
        self.columns = None
        self.mmi_data = None
        if output_dir is None:
            self.output_dir = os.path.dirname(grid_xml_path)
        else:
            self.output_dir = output_dir
        if output_basename is None:
            self.output_basename = "mmi"
        else:
            self.output_basename = output_basename
        self.algorithm_name = algorithm_filename_flag
        self.grid_xml_path = grid_xml_path
        self.parse_grid_xml()

    def extract_date_time(self, the_time_stamp):
        """Extract the parts of a date given a timestamp as per below example.

        :param the_time_stamp: The 'event_timestamp' attribute from  grid.xml.
        :type the_time_stamp: str

        # now separate out its parts
        # >>> e = "2012-08-07T01:55:12WIB"
        #>>> e[0:10]
        #'2012-08-07'
        #>>> e[12:-3]
        #'01:55:11'
        #>>> e[-3:]
        #'WIB'   (WIB = Western Indonesian Time)
        """
        date_tokens = the_time_stamp[0:10].split('-')
        self.year = int(date_tokens[0])
        self.month = int(date_tokens[1])
        self.day = int(date_tokens[2])

        time_tokens = the_time_stamp[11:-3].split(':')
        self.hour = int(time_tokens[0])
        self.minute = int(time_tokens[1])
        self.second = int(time_tokens[2])

        # right now only handles Indonesian Timezones
        tz_dict = {
            'WIB': timezone('Asia/Jakarta'),
            'WITA': timezone('Asia/Makassar'),
            'WIT': timezone('Asia/Jayapura')
        }
        tzinfo = tz_dict.get(self.time_zone)
        self.time = datetime(
            self.year,
            self.month,
            self.day,
            self.hour,
            self.minute,
            self.second,
            # For now realtime always uses Indonesia Time
            tzinfo=tzinfo
        )

    def parse_grid_xml(self):
        """Parse the grid xyz and calculate the bounding box of the event.

        :raises: GridXmlParseError

        The grid xyz dataset looks like this::

           <?xml version="1.0" encoding="US-ASCII" standalone="yes"?>
           <shakemap_grid xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
           xmlns="http://earthquake.usgs.gov/eqcenter/shakemap"
           xsi:schemaLocation="http://earthquake.usgs.gov
           http://earthquake.usgs.gov/eqcenter/shakemap/xml/schemas/
           shakemap.xsd" event_id="20120807015938" shakemap_id="20120807015938"
           shakemap_version="1" code_version="3.5"
           process_timestamp="2012-08-06T18:28:37Z" shakemap_originator="us"
           map_status="RELEASED" shakemap_event_type="ACTUAL">
           <event magnitude="5.1" depth="206" lat="2.800000"
               lon="128.290000" event_timestamp="2012-08-07T01:55:12WIB"
               event_network="" event_description="Halmahera, Indonesia    " />
           <grid_specification lon_min="126.290000" lat_min="0.802000"
               lon_max="130.290000" lat_max="4.798000"
               nominal_lon_spacing="0.025000" nominal_lat_spacing="0.024975"
               nlon="161" nlat="161" />
           <grid_field index="1" name="LON" units="dd" />
           <grid_field index="2" name="LAT" units="dd" />
           <grid_field index="3" name="PGA" units="pctg" />
           <grid_field index="4" name="PGV" units="cms" />
           <grid_field index="5" name="MMI" units="intensity" />
           <grid_field index="6" name="PSA03" units="pctg" />
           <grid_field index="7" name="PSA10" units="pctg" />
           <grid_field index="8" name="PSA30" units="pctg" />
           <grid_field index="9" name="STDPGA" units="pctg" />
           <grid_field index="10" name="URAT" units="" />
           <grid_field index="11" name="SVEL" units="ms" />
           <grid_data>
           126.2900 04.7980 0.01 0.02 1.16 0.05 0.02 0 0.5 1 600
           126.3150 04.7980 0.01 0.02 1.16 0.05 0.02 0 0.5 1 600
           126.3400 04.7980 0.01 0.02 1.17 0.05 0.02 0 0.5 1 600
           126.3650 04.7980 0.01 0.02 1.17 0.05 0.02 0 0.5 1 600
           ...
           ... etc

        .. note:: We could have also obtained some of this data from the
           grid.xyz and event.xml but the **grid.xml** is preferred because it
           contains clear and unequivical metadata describing the various
           fields and attributes. Also it provides all the data we need in a
           single file.
        """
        LOGGER.debug('ParseGridXml requested.')
        grid_path = self.grid_file_path()
        try:
            document = minidom.parse(grid_path)
            event_element = document.getElementsByTagName('event')
            event_element = event_element[0]
            self.magnitude = float(
                event_element.attributes['magnitude'].nodeValue)
            self.longitude = float(
                event_element.attributes['lon'].nodeValue)
            self.latitude = float(
                event_element.attributes['lat'].nodeValue)
            self.location = event_element.attributes[
                'event_description'].nodeValue.strip()
            self.depth = float(
                event_element.attributes['depth'].nodeValue)
            # Get the date - it's going to look something like this:
            # 2012-08-07T01:55:12WIB
            time_stamp = event_element.attributes['event_timestamp'].nodeValue
            # Note the timezone here is inconsistent with YZ from grid.xml
            # use the latter
            self.time_zone = time_stamp[19:]
            self.extract_date_time(time_stamp)

            specification_element = document.getElementsByTagName(
                'grid_specification')
            specification_element = specification_element[0]
            self.x_minimum = float(
                specification_element.attributes['lon_min'].nodeValue)
            self.x_maximum = float(
                specification_element.attributes['lon_max'].nodeValue)
            self.y_minimum = float(
                specification_element.attributes['lat_min'].nodeValue)
            self.y_maximum = float(
                specification_element.attributes['lat_max'].nodeValue)
            self.grid_bounding_box = QgsRectangle(
                self.x_minimum, self.y_maximum, self.x_maximum, self.y_minimum)
            self.rows = float(
                specification_element.attributes['nlat'].nodeValue)
            self.columns = float(
                specification_element.attributes['nlon'].nodeValue)
            data_element = document.getElementsByTagName(
                'grid_data')
            data_element = data_element[0]
            data = data_element.firstChild.nodeValue

            # Extract the 1,2 and 5th (MMI) columns and populate mmi_data
            longitude_column = 0
            latitude_column = 1
            mmi_column = 4
            self.mmi_data = []
            for line in data.split('\n'):
                if not line:
                    continue
                tokens = line.split(' ')
                longitude = tokens[longitude_column]
                latitude = tokens[latitude_column]
                mmi = tokens[mmi_column]
                mmi_tuple = (longitude, latitude, mmi)
                self.mmi_data.append(mmi_tuple)

        except Exception, e:
            LOGGER.exception('Event parse failed')
            raise GridXmlParseError(
                'Failed to parse grid file.\n%s\n%s' % (e.__class__, str(e)))

    def grid_file_path(self):
        """Validate that grid file path points to a file.

        :return: The grid xml file path.
        :rtype: str

        :raises: GridXmlFileNotFoundError
        """
        if os.path.isfile(self.grid_xml_path):
            return self.grid_xml_path
        else:
            raise GridXmlFileNotFoundError

    def mmi_to_delimited_text(self):
        """Return the mmi data as a delimited test string.

        :returns: A delimited text string that can easily be written to disk
            for e.g. use by gdal_grid.
        :rtype: str

        The returned string will look like this::

           123.0750,01.7900,1
           123.1000,01.7900,1.14
           123.1250,01.7900,1.15
           123.1500,01.7900,1.16
           etc...
        """
        delimited_text = 'lon,lat,mmi\n'
        for row in self.mmi_data:
            delimited_text += '%s,%s,%s\n' % (row[0], row[1], row[2])
        return delimited_text

    def mmi_to_delimited_file(self, force_flag=True):
        """Save mmi_data to delimited text file suitable for gdal_grid.

        The output file will be of the same format as strings returned from
        :func:`mmi_to_delimited_text`.

        :param force_flag: Whether to force the regeneration of the output
            file. Defaults to False.
        :type force_flag: bool

        :returns: The absolute file system path to the delimited text file.
        :rtype: str

        .. note:: An accompanying .csvt will be created which gdal uses to
           determine field types. The csvt will contain the following string:
           "Real","Real","Real". These types will be used in other conversion
           operations. For example to convert the csv to a shp you would do::

              ogr2ogr -select mmi -a_srs EPSG:4326 mmi.shp mmi.vrt mmi
        """
        LOGGER.debug('mmi_to_delimited_text requested.')

        csv_path = os.path.join(
            self.output_dir, 'mmi.csv')
        # short circuit if the csv is already created.
        if os.path.exists(csv_path) and force_flag is not True:
            return csv_path
        csv_file = file(csv_path, 'w')
        csv_file.write(self.mmi_to_delimited_text())
        csv_file.close()

        # Also write the .csvt which contains metadata about field types
        csvt_path = os.path.join(
            self.output_dir, self.output_basename + '.csvt')
        csvt_file = file(csvt_path, 'w')
        csvt_file.write('"Real","Real","Real"')
        csvt_file.close()

        return csv_path

    def mmi_to_vrt(self, force_flag=True):
        """Save the mmi_data to an ogr vrt text file.

        :param force_flag: Whether to force the regeneration of the output
            file. Defaults to False.
        :type force_flag: bool

        :returns: The absolute file system path to the .vrt text file.
        :rtype: str

        :raises: None
        """
        # Ensure the delimited mmi file exists
        LOGGER.debug('mmi_to_vrt requested.')

        vrt_path = os.path.join(
            self.output_dir,
            self.output_basename + '.vrt')

        # short circuit if the vrt is already created.
        if os.path.exists(vrt_path) and force_flag is not True:
            return vrt_path

        csv_path = self.mmi_to_delimited_file(True)

        vrt_string = (
            '<OGRVRTDataSource>'
            '  <OGRVRTLayer name="mmi">'
            '    <SrcDataSource>%s</SrcDataSource>'
            '    <GeometryType>wkbPoint</GeometryType>'
            '    <GeometryField encoding="PointFromColumns"'
            '                      x="lon" y="lat" z="mmi"/>'
            '  </OGRVRTLayer>'
            '</OGRVRTDataSource>' % csv_path)

        with codecs.open(vrt_path, 'w', encoding='utf-8') as f:
            f.write(vrt_string)

        return vrt_path

    def _add_executable_prefix(self, command):
        """Add the executable prefix for gdal binaries.

        This is primarily needed for OSX where gdal tools are tucked away in
        the Library path.

        :param command: A string containing the command to
        which the prefix will be prepended
        :type command: str

        :return: A copy of the command with the prefix added.
        :rtype: str

        :raises: None
        """

        executable_prefix = ''
        if sys.platform == 'darwin':  # Mac OS X
            # .. todo:: FIXME - softcode gdal version in this path
            executable_prefix = ('/Library/Frameworks/GDAL.framework/'
                                 'Versions/1.9/Programs/')
        command = executable_prefix + command
        return command

    def _run_command(self, command):
        """Run a command and raise any error as needed.

        This is a simple runner for executing gdal commands.

        :param command: A command string to be run.
        :type command: str

        :raises: Any exceptions will be propagated.
        """
        command = self._add_executable_prefix(command)
        try:
            my_result = call(command, shell=True)
            del my_result
        except CalledProcessError, e:
            LOGGER.exception('Running command failed %s' % command)
            message = (
                'Error while executing the following shell '
                'command: %s\nError message: %s' % (command, str(e)))
            # shameless hack - see https://github.com/AIFDR/inasafe/issues/141
            if sys.platform == 'darwin':  # Mac OS X
                if 'Errno 4' in str(e):
                    # continue as the error seems to be non critical
                    pass
                else:
                    raise Exception(message)
            else:
                raise Exception(message)

    def mmi_to_raster(
            self, force_flag=False, algorithm='nearest'):
        """Convert the grid.xml's mmi column to a raster using gdal_grid.

        A geotiff file will be created.

        Unfortunately no python bindings exist for doing this so we are
        going to do it using a shell call.

        .. see also:: http://www.gdal.org/gdal_grid.html

        Example of the gdal_grid call we generate::

           gdal_grid -zfield "mmi" -a invdist:power=2.0:smoothing=1.0 \
           -txe 126.29 130.29 -tye 0.802 4.798 -outsize 400 400 -of GTiff \
           -ot Float16 -l mmi mmi.vrt mmi.tif

        .. note:: It is assumed that gdal_grid is in your path.

        :param force_flag: Whether to force the regeneration of the output
            file. Defaults to False.
        :type force_flag: bool

        :param algorithm: Which re-sampling algorithm to use.
            valid options are 'nearest' (for nearest neighbour), 'invdist'
            (for inverse distance), 'average' (for moving average). Defaults
            to 'nearest' if not specified. Note that passing re-sampling alg
            parameters is currently not supported. If None is passed it will
            be replaced with 'nearest'.
        :type algorithm: str

        :returns: Path to the resulting tif file.
        :rtype: str

        .. note:: For interest you can also make quite beautiful smoothed
          raster using this:

          gdal_grid -zfield "mmi" -a_srs EPSG:4326
          -a invdist:power=2.0:smoothing=1.0 -txe 122.45 126.45
          -tye -2.21 1.79 -outsize 400 400 -of GTiff
          -ot Float16 -l mmi mmi.vrt mmi-trippy.tif
        """
        LOGGER.debug('mmi_to_raster requested.')

        if algorithm is None:
            algorithm = 'nearest'

        if self.algorithm_name:
            tif_path = os.path.join(
                self.output_dir, '%s-%s.tif' % (
                    self.output_basename, algorithm))
        else:
            tif_path = os.path.join(
                self.output_dir, '%s.tif' % self.output_basename)
        # short circuit if the tif is already created.
        if os.path.exists(tif_path) and force_flag is not True:
            return tif_path

        # Ensure the vrt mmi file exists (it will generate csv too if needed)
        vrt_path = self.mmi_to_vrt(force_flag)

        # now generate the tif using default nearest neighbour interpolation
        # options. This gives us the same output as the mi.grd generated by
        # the earthquake server.

        if 'invdist' in algorithm:
            algorithm = 'invdist:power=2.0:smoothing=1.0'

        # (Sunni): I'm not sure how this 'mmi' will work
        # (Tim): Its the mapping to which field in the CSV contains the data
        #    to be gridded.
        command = ((
            '%(gdal_grid)s -a %(alg)s -zfield "mmi" -txe %(xMin)s '
            '%(xMax)s -tye %(yMin)s %(yMax)s -outsize %(dimX)i '
            '%(dimY)i -of GTiff -ot Float16 -a_srs EPSG:4326 -l mmi '
            '"%(vrt)s" "%(tif)s"') % {
                'gdal_grid': which('gdal_grid')[0],
                'alg': algorithm,
                'xMin': self.x_minimum,
                'xMax': self.x_maximum,
                'yMin': self.y_minimum,
                'yMax': self.y_maximum,
                'dimX': self.columns,
                'dimY': self.rows,
                'vrt': vrt_path,
                'tif': tif_path
            })

        LOGGER.info('Created this gdal command:\n%s' % command)
        # Now run GDAL warp scottie...
        self._run_command(command)

        # We will use keywords file name with simple algorithm name since it
        # will raise an error in windows related to having double colon in path
        if 'invdist' in algorithm:
            algorithm = 'invdist'

        # copy the keywords file from fixtures for this layer
        self.create_keyword_file(algorithm)

        # Lastly copy over the standard qml (QGIS Style file) for the mmi.tif
        if self.algorithm_name:
            qml_path = os.path.join(
                self.output_dir, '%s-%s.qml' % (
                    self.output_basename, algorithm))
        else:
            qml_path = os.path.join(
                self.output_dir, '%s.qml' % self.output_basename)
        qml_source_path = os.path.join(data_dir(), 'mmi.qml')
        shutil.copyfile(qml_source_path, qml_path)
        return tif_path

    def mmi_to_shapefile(self, force_flag=False):
        """Convert grid.xml's mmi column to a vector shp file using ogr2ogr.

        An ESRI shape file will be created.

        :param force_flag: bool (Optional). Whether to force the regeneration
            of the output file. Defaults to False.

        :return: Path to the resulting tif file.
        :rtype: str

        Example of the ogr2ogr call we generate::

           ogr2ogr -select mmi -a_srs EPSG:4326 mmi.shp mmi.vrt mmi

        .. note:: It is assumed that ogr2ogr is in your path.
        """
        LOGGER.debug('mmi_to_shapefile requested.')

        shp_path = os.path.join(
            self.output_dir, '%s-points.shp' % self.output_basename)
        # Short circuit if the tif is already created.
        if os.path.exists(shp_path) and force_flag is not True:
            return shp_path

        # Ensure the vrt mmi file exists (it will generate csv too if needed)
        vrt_path = self.mmi_to_vrt(force_flag)

        # now generate the tif using default interpolation options

        command = (
            ('ogr2ogr -overwrite -select mmi -a_srs EPSG:4326 '
             '%(shp)s %(vrt)s mmi') % {'shp': shp_path, 'vrt': vrt_path})

        LOGGER.info('Created this gdal command:\n%s' % command)
        # Now run GDAL warp scottie...
        # noinspection PyProtectedMember
        self._run_command(command)

        # Lastly copy over the standard qml (QGIS Style file) for the mmi.tif
        qml_path = os.path.join(
            self.output_dir, '%s-points.qml' % self.output_basename)
        source_qml = os.path.join(data_dir(), 'mmi-shape.qml')
        shutil.copyfile(source_qml, qml_path)
        return shp_path

    def mmi_to_contours(self, force_flag=True, algorithm='nearest'):
        """Extract contours from the event's tif file.

        Contours are extracted at a 0.5 MMI interval. The resulting file will
        be saved in the extract directory. In the easiest use case you can

        :param force_flag:  (Optional). Whether to force the
         regeneration of contour product. Defaults to False.
        :type force_flag: bool

        :param algorithm: (Optional) Which interpolation algorithm to
                  use to create the underlying raster. Defaults to 'nearest'.
        :type algorithm: str
         **Only enforced if theForceFlag is true!**

        :returns: An absolute filesystem path pointing to the generated
            contour dataset.
        :exception: ContourCreationError

         simply do::

           shake_grid = ShakeGrid()
           contour_path = shake_grid.mmi_to_contours()

        which will return the contour dataset for the latest event on the
        ftp server.
        """
        LOGGER.debug('mmi_to_contours requested.')
        # TODO: Use sqlite rather?
        output_file_base = os.path.join(
            self.output_dir,
            '%s-contours-%s.' % (self.output_basename, algorithm))
        output_file = output_file_base + 'shp'
        if os.path.exists(output_file) and force_flag is not True:
            return output_file
        elif os.path.exists(output_file):
            try:
                os.remove(output_file_base + 'shp')
                os.remove(output_file_base + 'shx')
                os.remove(output_file_base + 'dbf')
                os.remove(output_file_base + 'prj')
            except OSError:
                LOGGER.exception(
                    'Old contour files not deleted'
                    ' - this may indicate a file permissions issue.')

        tif_path = self.mmi_to_raster(force_flag, algorithm)
        # Based largely on
        # http://svn.osgeo.org/gdal/trunk/autotest/alg/contour.py
        driver = ogr.GetDriverByName('ESRI Shapefile')
        ogr_dataset = driver.CreateDataSource(output_file)
        if ogr_dataset is None:
            # Probably the file existed and could not be overriden
            raise ContourCreationError(
                'Could not create datasource for:\n%s. Check that the file '
                'does not already exist and that you do not have file system '
                'permissions issues' % output_file)
        layer = ogr_dataset.CreateLayer('contour')
        field_definition = ogr.FieldDefn('ID', ogr.OFTInteger)
        layer.CreateField(field_definition)
        field_definition = ogr.FieldDefn('MMI', ogr.OFTReal)
        layer.CreateField(field_definition)
        # So we can fix the x pos to the same x coord as centroid of the
        # feature so labels line up nicely vertically
        field_definition = ogr.FieldDefn('X', ogr.OFTReal)
        layer.CreateField(field_definition)
        # So we can fix the y pos to the min y coord of the whole contour so
        # labels line up nicely vertically
        field_definition = ogr.FieldDefn('Y', ogr.OFTReal)
        layer.CreateField(field_definition)
        # So that we can set the html hex colour based on its MMI class
        field_definition = ogr.FieldDefn('RGB', ogr.OFTString)
        layer.CreateField(field_definition)
        # So that we can set the label in it roman numeral form
        field_definition = ogr.FieldDefn('ROMAN', ogr.OFTString)
        layer.CreateField(field_definition)
        # So that we can set the label horizontal alignment
        field_definition = ogr.FieldDefn('ALIGN', ogr.OFTString)
        layer.CreateField(field_definition)
        # So that we can set the label vertical alignment
        field_definition = ogr.FieldDefn('VALIGN', ogr.OFTString)
        layer.CreateField(field_definition)
        # So that we can set feature length to filter out small features
        field_definition = ogr.FieldDefn('LEN', ogr.OFTReal)
        layer.CreateField(field_definition)

        tif_dataset = gdal.Open(tif_path, GA_ReadOnly)
        # see http://gdal.org/java/org/gdal/gdal/gdal.html for these options
        band = 1
        contour_interval = 0.5
        contour_base = 0
        fixed_level_list = []
        use_no_data_flag = 0
        no_data_value = -9999
        id_field = 0  # first field defined above
        elevation_field = 1  # second (MMI) field defined above
        try:
            gdal.ContourGenerate(
                tif_dataset.GetRasterBand(band),
                contour_interval,
                contour_base,
                fixed_level_list,
                use_no_data_flag,
                no_data_value,
                layer,
                id_field,
                elevation_field)
        except Exception, e:
            LOGGER.exception('Contour creation failed')
            raise ContourCreationError(str(e))
        finally:
            del tif_dataset
            ogr_dataset.Release()

        # Copy over the standard .prj file since ContourGenerate does not
        # create a projection definition
        qml_path = os.path.join(
            self.output_dir,
            '%s-contours-%s.prj' % (self.output_basename, algorithm))
        source_qml = os.path.join(data_dir(), 'mmi-contours.prj')
        shutil.copyfile(source_qml, qml_path)

        # Lastly copy over the standard qml (QGIS Style file)
        qml_path = os.path.join(
            self.output_dir,
            '%s-contours-%s.qml' % (self.output_basename, algorithm))
        source_qml = os.path.join(data_dir(), 'mmi-contours.qml')
        shutil.copyfile(source_qml, qml_path)

        # Now update the additional columns - X,Y, ROMAN and RGB
        try:
            self.set_contour_properties(output_file)
        except InvalidLayerError:
            raise

        return output_file

    def set_contour_properties(self, input_file):
        """Set the X, Y, RGB, ROMAN attributes of the contour layer.

        :param input_file: (Required) Name of the contour layer.
        :type input_file: str

        :raise: InvalidLayerError if anything is amiss with the layer.
        """
        LOGGER.debug('set_contour_properties requested for %s.' % input_file)
        layer = QgsVectorLayer(input_file, 'mmi-contours', "ogr")
        if not layer.isValid():
            raise InvalidLayerError(input_file)

        layer.startEditing()
        # Now loop through the db adding selected features to mem layer
        request = QgsFeatureRequest()
        fields = layer.dataProvider().fields()

        for feature in layer.getFeatures(request):
            if not feature.isValid():
                LOGGER.debug('Skipping feature')
                continue
            # Work out x and y
            line = feature.geometry().asPolyline()
            y = line[0].y()

            x_max = line[0].x()
            x_min = x_max
            for point in line:
                if point.y() < y:
                    y = point.y()
                x = point.x()
                if x < x_min:
                    x_min = x
                if x > x_max:
                    x_max = x
            x = x_min + ((x_max - x_min) / 2)

            # Get length
            length = feature.geometry().length()

            mmi_value = float(feature['MMI'])
            # We only want labels on the whole number contours
            if mmi_value != round(mmi_value):
                roman = ''
            else:
                roman = romanise(mmi_value)

            # RGB from http://en.wikipedia.org/wiki/Mercalli_intensity_scale
            rgb = mmi_colour(mmi_value)

            # Now update the feature
            feature_id = feature.id()
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('X'), x)
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('Y'), y)
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('RGB'), rgb)
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('ROMAN'), roman)
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('ALIGN'), 'Center')
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('VALIGN'), 'HALF')
            layer.changeAttributeValue(
                feature_id, fields.indexFromName('LEN'), length)

        layer.commitChanges()

    def create_keyword_file(self, algorithm):
        """Create keyword file for the raster file created.

        Basically copy a template from keyword file in converter data
        and add extra keyword (usually a title)

        :param algorithm: Which re-sampling algorithm to use.
            valid options are 'nearest' (for nearest neighbour), 'invdist'
            (for inverse distance), 'average' (for moving average). Defaults
            to 'nearest' if not specified. Note that passing re-sampling alg
            parameters is currently not supported. If None is passed it will
            be replaced with 'nearest'.
        :type algorithm: str
        """
        if self.algorithm_name:
            keyword_path = os.path.join(
                self.output_dir, '%s-%s.keywords' % (
                    self.output_basename, algorithm))
        else:
            keyword_path = os.path.join(
                self.output_dir, '%s.keywords' % self.output_basename)
        mmi_keywords = os.path.join(data_dir(), 'mmi.keywords')
        shutil.copyfile(mmi_keywords, keyword_path)
        # append title and source to the keywords file
        if len(self.title.strip()) == 0:
            keyword_title = self.output_basename
        else:
            keyword_title = self.title
        with open(keyword_path, 'a') as keyword_file:
            keyword_file.write(get_string('title: %s \n' % keyword_title))
            keyword_file.write(get_string('source: %s ' % self.source))


def convert_mmi_data(
        grid_xml_path,
        title,
        source,
        output_path=None,
        algorithm=None,
        algorithm_filename_flag=True):
    """Convenience function to convert a single file.

    :param grid_xml_path: Path to the xml shake grid file.
    :type grid_xml_path: str

    :param title: The title of the earthquake.
    :type title: str

    :param source: The source of the shake data.
    :type source: str

    :param output_path: Specify which path to use as an alternative to the
        default.
    :type output_path: str

    :param algorithm: Type of algorithm to be used.
    :type algorithm: str

    :param algorithm_filename_flag: Flag whether to use the algorithm in the
        output file's name.
    :type algorithm_filename_flag: bool

    :returns: A path to the resulting raster file.
    :rtype: str
    """
    LOGGER.debug(grid_xml_path)
    LOGGER.debug(output_path)
    if output_path is not None:
        output_dir, output_basename = os.path.split(output_path)
        output_basename, _ = os.path.splitext(output_basename)
        LOGGER.debug('output_dir : ' + output_dir +
                     'output_basename : ' + output_basename)
    else:
        output_dir = output_path
        output_basename = None
    converter = ShakeGrid(
        title=title,
        source=source,
        grid_xml_path=grid_xml_path,
        output_dir=output_dir,
        output_basename=output_basename,
        algorithm_filename_flag=algorithm_filename_flag)
    return converter.mmi_to_raster(force_flag=True, algorithm=algorithm)
