# coding=utf-8
"""**Convert QgsRasterLayer to QgsVectorLayer**
"""
__author__ = 'Dmitry Kolesov <kolesov.dm@gmail.com>'
__revision__ = '502f00bf1e6755381fba9bb8b68ee10c70035e70'
__date__ = '14/01/2014'
__license__ = "GPL"
__copyright__ = 'Copyright 2012, Australia Indonesia Facility for '
__copyright__ += 'Disaster Reduction'


from qgis.core import (
    QgsRasterLayer,
    QgsRectangle,
    QgsField,
    QgsVectorLayer,
    QgsFeature,
    QgsPoint,
    QgsGeometry,
    QgsRasterFileWriter,
    QgsRasterPipe
)
# noinspection PyPackageRequirements
from PyQt4.QtCore import QVariant
from qgis_vector_tools import (
    union_geometry,
    points_to_rectangles
)

from safe.common.utilities import unique_filename
from safe.gis.gdal_ogr_tools import polygonize_thresholds
from safe.common.exceptions import GetDataError


def _get_pixel_coordinates(extent, width, height, row, col):
    """Pixel to coordinates transformation

    :param extent: Geotransformation parameter: extent
    :type extent: QgsRectangle

    :param width: Geotransformation parameter: horizontal raster size
    :type width: Int

    :param height: Geotransformation parameter: vertical raster size
    :type height: Int

    :param col: Input x pixel coordinate
    :type col: Int

    :param row: Input y pixel coordinate
    :type row: Int

    :returns: output_x,output_y Output coordinates
    :rtype: (double, double)
    """
    # pylint: disable=
    x_minimum = extent.xMinimum()
    y_maximum = extent.yMaximum()
    x_resolution = (extent.xMaximum() - extent.xMinimum()) / width
    y_resolution = (extent.yMaximum() - extent.yMinimum()) / height

    output_x = x_minimum + col * x_resolution
    output_y = y_maximum - row * y_resolution
    return output_x, output_y


def pixels_to_points(
        raster,
        threshold_min=0.0,
        threshold_max=float('inf'),
        field_name='value'):
    """
    Convert raster to points.

    Areas (pixels) with threshold_min < pixel_values < threshold_max will be
    converted to point layer.

    :param raster: Raster layer
    :type raster: QgsRasterLayer

    :param threshold_min: Value that splits raster to flooded or not flooded.
    :type threshold_min: float

    :param threshold_max: Value that splits raster to flooded or not flooded.
    :type threshold_max: float

    :param field_name: Field name to store pixel value.
    :type field_name:  string

    :returns: Point layer of pixels
    :rtype: QgsVectorLayer
    """
    if raster.bandCount() != 1:
        msg = "Current version allows using of one-band raster only"
        raise NotImplementedError(msg)

    extent = raster.extent()
    width, height = raster.width(), raster.height()
    provider = raster.dataProvider()
    block = provider.block(1, extent, width, height)

    # Create points
    crs = raster.crs().toWkt()
    point_layer = QgsVectorLayer('Point?crs=' + crs, 'pixels', 'memory')

    point_provider = point_layer.dataProvider()
    point_provider.addAttributes([QgsField(field_name, QVariant.Double)])
    field_index = point_provider.fieldNameIndex(field_name)
    attribute_count = 1

    point_layer.startEditing()
    for row in range(height):
        for col in range(width):
            value = block.value(row, col)
            x, y = _get_pixel_coordinates(extent, width, height, row, col)
            # noinspection PyCallByClass,PyTypeChecker,PyArgumentList
            geom = QgsGeometry.fromPoint(QgsPoint(x, y))
            if threshold_min < value < threshold_max:
                feature = QgsFeature()
                feature.initAttributes(attribute_count)
                feature.setAttribute(field_index, value)
                feature.setGeometry(geom)
                _ = point_layer.dataProvider().addFeatures([feature])
    point_layer.commitChanges()
    return point_layer


def polygonize(raster, threshold_min=0.0, threshold_max=float('inf')):
    """Raster polygonizer.

    Function to polygonize raster. Areas (pixels) with threshold_min <
    pixel_values < threshold_max will be converted to polygons.

    :param raster: Raster layer
    :type raster: QgsRasterLayer

    :param threshold_min: Value that splits raster to flooded or not flooded.
    :type threshold_min: float

    :param threshold_max: Value that splits raster to flooded or not flooded.
    :type threshold_max: float

    :returns: Polygonal geometry
    :rtype: QgsGeometry
    """
    points = pixels_to_points(raster, threshold_min, threshold_max)
    polygons = points_to_rectangles(
        points,
        raster.rasterUnitsPerPixelX(),
        raster.rasterUnitsPerPixelY())
    polygons = union_geometry(polygons)
    return polygons


def align_clip_raster(hazard_layer, extent):
    """Align raster extent and the given extent.

    Both layers should be in the same projection.

    .. versionadded:: 3.4
    .. note:: Delegates to clip_raster()

    :param hazard_layer: The raster layer.
    :type hazard_layer: QgsRasterLayer

    :param extent: The extent in the form [xmin, ymin, xmax, ymax]
    :type extent: list

    :returns: Clipped region of the raster
    :rtype: QgsRasterLayer
    """
    raster_extent = hazard_layer.dataProvider().extent()
    clip_xmin = raster_extent.xMinimum()
    clip_ymin = raster_extent.yMinimum()

    if extent[0] > clip_xmin:
        clip_xmin = extent[0]
    if extent[1] > clip_ymin:
        clip_ymin = extent[1]
    height = ((extent[3] - extent[1]) /
              hazard_layer.rasterUnitsPerPixelY())
    height = int(height)
    width = ((extent[2] - extent[0]) /
             hazard_layer.rasterUnitsPerPixelX())
    width = int(width)
    raster_extent = hazard_layer.dataProvider().extent()
    xmin = raster_extent.xMinimum()
    xmax = raster_extent.xMaximum()
    ymin = raster_extent.yMinimum()
    ymax = raster_extent.yMaximum()
    x_delta = (xmax - xmin) / hazard_layer.width()
    x = xmin
    for i in range(hazard_layer.width()):
        if abs(x - clip_xmin) < x_delta:
            # We have found the aligned raster boundary
            break
        x += x_delta
        _ = i
    y_delta = (ymax - ymin) / hazard_layer.height()
    y = ymin
    for i in range(hazard_layer.width()):
        if abs(y - clip_ymin) < y_delta:
            # We have found the aligned raster boundary
            break
        y += y_delta
    clip_extent = [x, y, x + width * x_delta, y + height * y_delta]
    # Clip hazard raster
    small_raster = clip_raster(
        hazard_layer, width, height, QgsRectangle(*clip_extent))
    return small_raster


def clip_raster(raster, column_count, row_count, output_extent):
    """Clip raster to specified extent, width and height.

    Note there is similar utility in safe_qgis.utilities.clipper, but it uses
    gdal whereas this one uses native QGIS.

    :param raster: Raster
    :type raster: QgsRasterLayer

    :param column_count: Desired width in pixels of new raster
    :type column_count: Int

    :param row_count: Desired height in pixels of new raster
    :type row_count: Int

    :param output_extent: Extent of the clipped region
    :type output_extent: QgsRectangle

    :returns: Clipped region of the raster
    :rtype: QgsRasterLayer
    """
    provider = raster.dataProvider()
    pipe = QgsRasterPipe()
    pipe.set(provider.clone())

    base_name = unique_filename()
    file_name = base_name + '.tif'
    file_writer = QgsRasterFileWriter(file_name)
    file_writer.writeRaster(
        pipe,
        column_count,
        row_count,
        output_extent,
        raster.crs())

    return QgsRasterLayer(file_name, 'clipped_raster')


def polygonize_gdal(
        raster,
        threshold_min=0.0,
        threshold_max=float('inf')):
    """
    Function to polygonize raster. Areas (pixels) with threshold_min <
    pixel_values < threshold_max will be converted to polygons.

    :param raster:  Raster layer
    :type raster: QgsRasterLayer

    :param threshold_min: Value that splits raster to flooded or not flooded.
    :type threshold_min: float

    :param threshold_max: Value that splits raster to flooded or not flooded.
    :type threshold_max: float

    :returns: Polygonal geometry
    :rtype: QgsGeometry
    """

    # save qgis raster to disk
    base_name = unique_filename()
    file_name = base_name + '.tif'
    file_writer = QgsRasterFileWriter(file_name)
    pipe = QgsRasterPipe()
    provider = raster.dataProvider()
    if not pipe.set(provider.clone()):
        msg = "Cannot set pipe provider"
        raise GetDataError(msg)
    file_writer.writeRaster(
        pipe,
        provider.xSize(),
        provider.ySize(),
        provider.extent(),
        provider.crs())

    (
        inside_file_name,
        inside_layer_name,
        outside_file_name,
        outside_layer_name
    ) = polygonize_thresholds(file_name, threshold_min, threshold_max)
    inside_layer = \
        QgsVectorLayer(inside_file_name, inside_layer_name, 'ogr')
    outside_layer = \
        QgsVectorLayer(outside_file_name, outside_layer_name, 'ogr')
    if inside_layer.featureCount() == 0:
        return None, None
    else:
        return inside_layer, outside_layer
