# coding=utf-8
"""A converter for USGS shakemap grid.xml files."""


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

import numpy as np
import pytz
# This import is required to enable PyQt API v2
# noinspection PyUnresolvedReferences
import qgis  # NOQA pylint: disable=unused-import
from osgeo import gdal, ogr
from osgeo.gdalconst import GA_ReadOnly
from pytz import timezone
from qgis.core import (
    QgsRectangle,
    QgsRasterLayer)

from safe.common.exceptions import (
    GridXmlFileNotFoundError,
    GridXmlParseError,
    ContourCreationError,
    InvalidLayerError,
    CallGDALError)
from safe.common.utilities import which
from safe.definitions.constants import (
    NONE_SMOOTHING, NUMPY_SMOOTHING, SCIPY_SMOOTHING
)
from safe.definitions.exposure import exposure_all
from safe.definitions.extra_keywords import (
    extra_keyword_earthquake_latitude,
    extra_keyword_earthquake_longitude,
    extra_keyword_earthquake_magnitude,
    extra_keyword_earthquake_depth,
    extra_keyword_earthquake_description,
    extra_keyword_earthquake_location,
    extra_keyword_time_zone,
    extra_keyword_earthquake_x_maximum,
    extra_keyword_earthquake_x_minimum,
    extra_keyword_earthquake_y_maximum,
    extra_keyword_earthquake_y_minimum,
    extra_keyword_earthquake_event_time,
    extra_keyword_earthquake_event_id,
)
from safe.definitions.hazard import hazard_earthquake
from safe.definitions.hazard_category import hazard_category_single_event
from safe.definitions.hazard_classifications import (
    earthquake_mmi_scale, generic_hazard_classes)
from safe.definitions.layer_geometry import layer_geometry_raster
from safe.definitions.layer_modes import layer_mode_continuous
from safe.definitions.layer_purposes import layer_purpose_hazard
from safe.definitions.units import unit_mmi
from safe.definitions.utilities import default_classification_thresholds
from safe.definitions.versions import inasafe_keyword_version
from safe.gis.raster.contour import (
    gaussian_kernel, convolve, set_contour_properties)
from safe.utilities.i18n import tr
from safe.utilities.keyword_io import KeywordIO
from safe.utilities.resources import resources_path

__copyright__ = "Copyright 2017, The InaSAFE Project"
__license__ = "GPL version 3"
__email__ = "info@inasafe.org"
__revision__ = '$Format:%H$'

LOGGER = logging.getLogger('InaSAFE')

# Algorithms
USE_ASCII = 'use_ascii'
NEAREST_NEIGHBOUR = 'nearest'
INVDIST = 'invdist'


class ShakeGrid():

    """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,
            smoothing_method=NONE_SMOOTHING,
            smooth_sigma=0.9,
            extra_keywords=None
    ):
        """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

        :param smoothing_method: Smoothing method. One of None, Numpy, Scipy.
        :type smoothing_method: basestring

        :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
        self.event_id = 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.smoothing_method = smoothing_method
        self.smoothing_sigma = smooth_sigma
        self.extra_keywords = extra_keywords if extra_keywords else {}

        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:19].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': 'Asia/Jakarta',
            'WITA': 'Asia/Makassar',
            'WIT': 'Asia/Jayapura'
        }
        if self.time_zone in tz_dict:
            self.time_zone = tz_dict.get(self.time_zone, self.time_zone)

        # noinspection PyBroadException
        try:
            if not self.time_zone:
                # default to utc if empty
                tzinfo = pytz.utc
            else:
                tzinfo = timezone(self.time_zone)
        except BaseException:
            tzinfo = pytz.utc

        self.time = datetime(
            self.year,
            self.month,
            self.day,
            self.hour,
            self.minute,
            self.second)
        # For now realtime always uses Western Indonesia Time
        self.time = tzinfo.localize(self.time)

    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)
            shakemap_grid_element = document.getElementsByTagName(
                'shakemap_grid')
            shakemap_grid_element = shakemap_grid_element[0]
            self.event_id = shakemap_grid_element.attributes[
                'event_id'].nodeValue

            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 = int(float(
                specification_element.attributes['nlat'].nodeValue))
            self.columns = int(float(
                specification_element.attributes['nlon'].nodeValue))
            data_element = document.getElementsByTagName(
                'grid_data')
            data_element = data_element[0]
            data = data_element.firstChild.nodeValue

            # Find grid_field of LON LAT and MMI
            grid_field = ['LON', 'LAT', 'MMI']
            grid_fields = document.getElementsByTagName('grid_field')
            grid_fields = {
                f.attributes['name'].nodeValue: int(
                    f.attributes['index'].nodeValue)
                for f in grid_fields
                if f.attributes['name'].nodeValue in grid_field
            }
            # grid_field index are one-based
            longitude_column = grid_fields['LON'] - 1
            latitude_column = grid_fields['LAT'] - 1
            mmi_column = grid_fields['MMI'] - 1
            lon_list = []
            lat_list = []
            mmi_list = []
            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]
                lon_list.append(float(longitude))
                lat_list.append(float(latitude))
                mmi_list.append(float(mmi))

            if self.smoothing_method == NUMPY_SMOOTHING:
                LOGGER.debug('We are using NUMPY smoothing')
                ncols = len(np.where(np.array(lon_list) == lon_list[0])[0])
                nrows = len(np.where(np.array(lat_list) == lat_list[0])[0])

                # reshape mmi_list to 2D array to apply gaussian filter
                Z = np.reshape(mmi_list, (nrows, ncols))

                # smooth MMI matrix
                mmi_list = convolve(Z, gaussian_kernel(self.smoothing_sigma))

                # reshape array back to 1D long list of mmi
                mmi_list = np.reshape(mmi_list, ncols * nrows)

            elif self.smoothing_method == SCIPY_SMOOTHING:
                LOGGER.debug('We are using SCIPY smoothing')
                from scipy.ndimage.filters import gaussian_filter
                ncols = len(np.where(np.array(lon_list) == lon_list[0])[0])
                nrows = len(np.where(np.array(lat_list) == lat_list[0])[0])

                # reshape mmi_list to 2D array to apply gaussian filter
                Z = np.reshape(mmi_list, (nrows, ncols))

                # smooth MMI matrix
                # Help from Hadi Ghasemi
                mmi_list = gaussian_filter(Z, self.smoothing_sigma)

                # reshape array back to 1D long list of mmi
                mmi_list = np.reshape(mmi_list, ncols * nrows)

            # zip lists as list of tuples
            self.mmi_data = list(zip(lon_list, lat_list, mmi_list))

        except Exception as 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 = open(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 = open(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 _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.
        """
        try:
            my_result = call(command, shell=True)
            del my_result
        except CalledProcessError as 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=USE_ASCII):
        """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 'use_ascii'.
            'use_ascii' algorithm will convert the mmi grid to ascii file
            then convert it to raster using gdal_translate.
        :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 = USE_ASCII

        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

        if algorithm == USE_ASCII:
            # Convert to ascii
            ascii_path = self.mmi_to_ascii(True)

            # Creating command to convert to tif
            command = (
                (
                    '%(gdal_translate)s -a_srs EPSG:4326 '
                    '"%(ascii)s" "%(tif)s"'
                ) % {
                    'gdal_translate': which('gdal_translate')[0],
                    'ascii': ascii_path,
                    'tif': tif_path
                }
            )

            LOGGER.info('Created this gdal command:\n%s' % command)
            # Now run GDAL warp scottie...
            self._run_command(command)
        else:
            # 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
            # mmi.grd generated by the earthquake server.

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

            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 = resources_path('converter_data', '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

        binary_list = which('ogr2ogr')
        LOGGER.debug('Path for ogr2ogr: %s' % binary_list)
        if len(binary_list) < 1:
            raise CallGDALError(
                tr('ogr2ogr could not be found on your computer'))
        # Use the first matching gdalwarp found
        binary = binary_list[0]
        command = (
            ('%(ogr2ogr)s -overwrite -select mmi -a_srs EPSG:4326 '
             '%(shp)s %(vrt)s mmi') % {
                'ogr2ogr': binary,
                'shp': shp_path,
                'vrt': vrt_path})

        LOGGER.info('Created this ogr command:\n%s' % command)
        # Now run ogr2ogr ...
        # 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 = resources_path('converter_data', 'mmi-shape.qml')
        shutil.copyfile(source_qml, qml_path)
        return shp_path

    def mmi_to_contours(self, force_flag=True, algorithm=USE_ASCII):
        """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))
        # There are minor issues with shapefile, so we switch to gpkg
        # See https://github.com/inasafe/inasafe/issues/5063
        output_file = output_file_base + 'gpkg'
        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 + 'gpkg')
                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
        # Use Geopackage driver to overcome this:
        # See https://github.com/inasafe/inasafe/issues/5063
        driver = ogr.GetDriverByName('GPKG')
        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 as 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
        projection_path = os.path.join(
            self.output_dir,
            '%s-contours-%s.prj' % (self.output_basename, algorithm))
        source_projection_path = resources_path(
            'converter_data', 'mmi-contours.prj')
        shutil.copyfile(source_projection_path, projection_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_path = resources_path('converter_data', 'mmi-contours.qml')
        shutil.copyfile(source_qml_path, qml_path)

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

        return output_file

    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
        """
        keyword_io = KeywordIO()

        # Set thresholds for each exposure
        mmi_default_classes = default_classification_thresholds(
            earthquake_mmi_scale
        )
        mmi_default_threshold = {
            earthquake_mmi_scale['key']: {
                'active': True,
                'classes': mmi_default_classes
            }
        }
        generic_default_classes = default_classification_thresholds(
            generic_hazard_classes
        )
        generic_default_threshold = {
            generic_hazard_classes['key']: {
                'active': True,
                'classes': generic_default_classes
            }
        }

        threshold_keyword = {}
        for exposure in exposure_all:
            # Not all exposure is supported by earthquake_mmi_scale
            if exposure in earthquake_mmi_scale['exposures']:
                threshold_keyword[exposure['key']] = mmi_default_threshold
            else:
                threshold_keyword[
                    exposure['key']] = generic_default_threshold

        extra_keywords = {
            extra_keyword_earthquake_latitude['key']: self.latitude,
            extra_keyword_earthquake_longitude['key']: self.longitude,
            extra_keyword_earthquake_magnitude['key']: self.magnitude,
            extra_keyword_earthquake_depth['key']: self.depth,
            extra_keyword_earthquake_description['key']: self.description,
            extra_keyword_earthquake_location['key']: self.location,
            extra_keyword_earthquake_event_time['key']: self.time.strftime(
                '%Y-%m-%dT%H:%M:%S'),
            extra_keyword_time_zone['key']: self.time_zone,
            extra_keyword_earthquake_x_minimum['key']: self.x_minimum,
            extra_keyword_earthquake_x_maximum['key']: self.x_maximum,
            extra_keyword_earthquake_y_minimum['key']: self.y_minimum,
            extra_keyword_earthquake_y_maximum['key']: self.y_maximum,
            extra_keyword_earthquake_event_id['key']: self.event_id
        }
        for key, value in list(self.extra_keywords.items()):
            extra_keywords[key] = value

        # Delete empty element.
        empty_keys = []
        for key, value in list(extra_keywords.items()):
            if value is None:
                empty_keys.append(key)
        for empty_key in empty_keys:
            extra_keywords.pop(empty_key)
        keywords = {
            'hazard': hazard_earthquake['key'],
            'hazard_category': hazard_category_single_event['key'],
            'keyword_version': inasafe_keyword_version,
            'layer_geometry': layer_geometry_raster['key'],
            'layer_mode': layer_mode_continuous['key'],
            'layer_purpose': layer_purpose_hazard['key'],
            'continuous_hazard_unit': unit_mmi['key'],
            'classification': earthquake_mmi_scale['key'],
            'thresholds': threshold_keyword,
            'extra_keywords': extra_keywords,
            'active_band': 1
        }

        if self.algorithm_name:
            layer_path = os.path.join(
                self.output_dir, '%s-%s.tif' % (
                    self.output_basename, algorithm))
        else:
            layer_path = os.path.join(
                self.output_dir, '%s.tif' % self.output_basename)

        # append title and source to the keywords file
        if len(self.title.strip()) == 0:
            keyword_title = self.output_basename
        else:
            keyword_title = self.title

        keywords['title'] = keyword_title

        hazard_layer = QgsRasterLayer(layer_path, keyword_title)

        if not hazard_layer.isValid():
            raise InvalidLayerError()

        keyword_io.write_keywords(hazard_layer, keywords)

    def mmi_to_ascii(self, force_flag=False):
        """Convert grid.xml mmi column to a ascii raster file."""
        ascii_path = os.path.join(
            self.output_dir, '%s.asc' % self.output_basename)
        # Short circuit if the tif is already created.
        if os.path.exists(ascii_path) and force_flag is not True:
            return ascii_path

        cell_size = (self.x_maximum - self.x_minimum) / (self.rows - 1)
        asc_file = open(ascii_path, 'w')
        asc_file.write('ncols %d\n' % self.columns)
        asc_file.write('nrows %d\n' % self.rows)
        asc_file.write('xllcorner %.3f\n' % self.x_minimum)
        asc_file.write('yllcorner %.3f\n' % self.y_minimum)
        asc_file.write('cellsize %.3f\n' % cell_size)
        asc_file.write('nodata_value -9999\n')

        cell_string = ''
        cell_values = np.reshape(
            [v[2] for v in self.mmi_data], (self.rows, self.columns))
        for i in range(self.rows):
            for j in range(self.columns):
                cell_string += '%.3f ' % cell_values[i][j]
            cell_string += '\n'

        asc_file.write(cell_string)

        asc_file.close()

        return ascii_path


def convert_mmi_data(
        grid_xml_path,
        title,
        source,
        output_path=None,
        algorithm=None,
        algorithm_filename_flag=True,
        smoothing_method=NONE_SMOOTHING,
        smooth_sigma=0.9,
        extra_keywords=None
):
    """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 place_layer: Nearby cities/places.
    :type place_layer: QgsVectorLayer

    :param place_name: Column name that indicates name of the cities.
    :type place_name: str

    :param place_population: Column name that indicates number of population.
    :type place_population: 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

    :param smoothing_method: Smoothing method. One of None, Numpy, Scipy.
    :type smoothing_method: basestring

    :param smooth_sigma: parameter for gaussian filter used in smoothing
        function.
    :type smooth_sigma: float

    :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,
        source,
        grid_xml_path,
        output_dir=output_dir,
        output_basename=output_basename,
        algorithm_filename_flag=algorithm_filename_flag,
        smoothing_method=smoothing_method,
        smooth_sigma=smooth_sigma,
        extra_keywords=extra_keywords
    )
    return converter.mmi_to_raster(force_flag=True, algorithm=algorithm)
