import json
import itertools
import rtree
import requests
from parse import search
from geomet import wkt
from datetime import datetime, timedelta

from functools import partial, lru_cache
import pyproj
import numpy as np
from time import time
import logging

from shapely.geometry import shape, mapping, box
from shapely.geometry.polygon import Polygon
from shapely.geometry.multipolygon import MultiPolygon
from shapely.geometry.collection import GeometryCollection
from shapely.geometry.point import Point
from shapely.ops import unary_union, transform


__all__ = ['json2ogr', 'ogr2json', 'dissolve', 'intersect', 'project_local',
           'project_global', 'buffer_to_dist', 'get_area', 'get_area_percent',
           'esri_server2ogr', 'get_species_count', 'esri_server2histo',
           'esri_count_groupby', 'cartodb2ogr', 'esri_count_30days',
           'esri_last_instance', 'erase', 'get_date_from_timestamp',
           'get_feature_count', 'test_ip', 'esri_attributes', 'get_presence',
           'get_histo_loss_area', 'get_histo_pre2001_area', 'get_histo_total_area',
           'get_area_by_attributes', 'get_geom_by_attributes', 'pad_counts',
           'vals_by_year', 'split', 'split_features', 'get_counts_by_year',
           'get_count_by_year', 'combine_counts_by_year', 'get_ok']


HA_CONVERSION = 10000
FEATURE_GEOG_THRESHOLD = 1.2
REQUEST_GEOG_THRESHOLD = 20
REQUEST_COORD_THRESHOLD = 20000
FUNCTION_COUNT = 0


def get_ok():
    return 'ok'


def test_ip():
    return requests.get('http://checkip.amazonaws.com').text.replace('\n', '')


def json2ogr(in_json):
    '''
    Convert geojson object to GDAL geometry
    '''
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION json2ogr STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    if isinstance(in_json, str):
        in_json = json.loads(in_json)

    if not isinstance(in_json, dict):
        raise ValueError('input json must be dictionary')

    if 'features' not in in_json.keys():
        raise ValueError('input json must contain features property')

    for f in in_json['features']:
        f['geometry'] = shape(f['geometry'])
        if not f['geometry'].is_valid:
            f['geometry'] = f['geometry'].buffer(0)

    for i in range(len(in_json['features'])):
        in_json['features'][i]['properties']['id'] = i

    # logging.info('FUNCTION json2ogr STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return in_json


def ogr2json(featureset):
    '''
    Convert GDAL geometry to geojson object
    '''
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION ogr2json STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    new_features = []
    for f in featureset['features']:
        new_features.append(dict(geometry=mapping(f['geometry']),
                                 properties=f['properties'],
                                 type=f['type']))
        # f['geometry'] = mapping(f['geometry'])

    new_featureset = dict(type=featureset['type'],
                          features=new_features)
    if 'crs' in featureset.keys():
        new_featureset['crs'] = featureset['crs']
    # logging.info('FUNCTION ogr2json STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return json.dumps(new_featureset)


def explode(coords):
    """Explode a GeoJSON geometry's coordinates object and yield coordinate
    tuples. As long as the input is conforming, the type of the geometry
    doesn't matter.
    https://gis.stackexchange.com/questions/90553/fiona-get-each-feature-
    extent-bounds"""
    for e in coords:
        if isinstance(e, (float, int)):
            yield coords
            break
        else:
            for f in explode(e):
                yield f


def bounds(f):
    if isinstance(f['geometry'], dict):
        geom = f['geometry']['coordinates']
    else:
        try:
            geom = mapping(f['geometry'])['coordinates']
        except Exception as e:
            raise ValueError((str(e),f['geometry'],mapping(f['geometry'])))
    x, y = zip(*list(explode(geom)))
    return min(x), min(y), max(x), max(y)


def bbox(f):
    tups = mapping(box(*bounds(f)))['coordinates']
    return [[list(tup) for tup in tups[0]]]


def ogr2rings(f):
    return [[list(tup) for tup in mapping(f['geometry'])['coordinates'][0]]]


# @lru_cache(5)
def esri_server2ogr(layer_endpoint, aoi, out_fields, where='1=1', token=''):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION esri_server2ogr STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    url = layer_endpoint.replace('?f=pjson', '') + '/query'

    params = {}
    params['where'] = where
    if 'objectid' not in out_fields:
        out_fields = 'objectid,' + out_fields if out_fields else 'objectid'
    params['outFields'] = out_fields
    params['returnGeometry'] = True
    params['returnM'] = False
    params['returnZ'] = False
    params['f'] = 'geojson'
    params['geometryType'] = 'esriGeometryPolygon'
    params['spatialRel'] = 'esriSpatialRelIntersects'

    # if protected service, retrieve token
    if token:
        params['token'] = token

    # iterate through aoi features (Esri does not accept multipart polygons
    # as a spatial filter, and the aoi features may be too far apart to combine
    # into one bounding box)
    featureset = {'type': 'FeatureCollection', 'features': []}
    features = []
    objectids = []
    if isinstance(aoi, str):
        aoi = json.loads(aoi)
    for f in aoi['features']:
        params['geometry'] = str({'rings': bbox(f),
                                  'spatialReference': {'wkid': 4326}})
        req = requests.post(url, data=params)
        req.raise_for_status()
        try:
            # response = json2ogr(req.text)
            response = req.json()
            assert 'features' in response
        except:
            raise ValueError((req.text, url, params))

        # append response to full dataset, except features already included
        for h in response['features']:
            feat_id = ','.join([str(prop) for prop in h['properties'].values()])
            if feat_id not in objectids:
                features.append(h)
                objectids.append(feat_id)

    featureset = json2ogr(dict(type='FeatureCollection',
                               features=features))

    # logging.info('FUNCTION esri_server2ogr STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return featureset


def esri_server2histo(layer_endpoint, aoi):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION esri_server2histo STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    url = layer_endpoint.replace('?f=pjson', '') + '/computeHistograms'

    params = {}
    params['f'] = 'json'
    params['geometryType'] = 'esriGeometryPolygon'
    params['spatialRel'] = 'esriSpatialRelIntersects'
    params['returnGeometry'] = True
    params['where'] = '1=1'

    if isinstance(aoi, str):
        aoi = json.loads(aoi)
    # if featureset['features']:
    #     f = featureset['features'][0]
    histogram = [0] * 256
    for f in aoi['features']:
        params['geometry'] = str({'rings': ogr2rings(f),
                                  'spatialReference': {'wkid': 4326}})
        req = requests.post(url, data=params)
        req.raise_for_status()
        try:
            response = req.json()['histograms']
            if response:
                for i, count in enumerate(response[0]['counts']):
                    histogram[i] += count
        except Exception as e:
            raise ValueError('{} --- {}'.format(e, req.text))

    # logging.info('FUNCTION esri_server2histo STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return histogram


def esri_attributes(layer_endpoint, aoi, out_fields):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION esri_attributes STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    url = layer_endpoint.replace('?f=pjson', '') + '/query'

    params = {}
    params['f'] = 'json'
    params['geometryType'] = 'esriGeometryPolygon'
    params['where'] = '1=1'
    params['spatialRel'] = 'esriSpatialRelIntersects'
    params['returnGeometry'] = False
    params['outFields'] = out_fields

    if isinstance(aoi, str):
        aoi = json.loads(aoi)
    objectids = []
    attributes = []
    for f in aoi['features']:
        params['geometry'] = str({'rings': ogr2rings(f),
                                  'spatialReference': {'wkid': 4326}})
        req = requests.post(url, data=params)
        req.raise_for_status()

        # return [feat['attributes'] for feat in req.json()['features']]
        for h in req.json()['features']:
            feat_id = ','.join([str(prop) for prop in h['attributes'].values()])
            if feat_id not in objectids:
                attributes.append(h['attributes'])
                objectids.append(feat_id)

    # logging.info('FUNCTION esri_attributes STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return attributes


def esri_count_groupby(layer_endpoint, aoi, count_fields):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION esri_count_groupby STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    url = layer_endpoint.replace('?f=pjson', '') + '/query'

    params = {}
    params['f'] = 'json'
    params['geometryType'] = 'esriGeometryPolygon'
    params['where'] = '1=1'
    params['spatialRel'] = 'esriSpatialRelIntersects'
    params['returnGeometry'] = False
    params['groupByFieldsForStatistics'] = count_fields
    count_fields = count_fields.split(',')
    params['outStatistics'] = json.dumps([{
        'statisticType': 'count',
        'onStatisticField': count_fields[0],
        'outStatisticFieldName': 'count'
    }])

    if isinstance(aoi, str):
        aoi = json.loads(aoi)
    counts = {}
    for f in aoi['features']:
        params['geometry'] = str({'rings': ogr2rings(f),
                                  'spatialReference': {'wkid': 4326}})
        req = requests.post(url, data=params)
        req.raise_for_status()
        try:
            f_counts = {'-'.join([str(item['attributes'][field]) for field in
                                  count_fields]): item['attributes']['count']
                        for item in req.json()['features']}
            for key, val in f_counts.items():
                if not key in counts.keys():
                    counts[key] = val
                else:
                    counts[key] += val

        except Exception as e:
            raise ValueError((str(e), url, params, req.text))

    # logging.info('FUNCTION esri_count_groupby STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return counts


def esri_count_30days(layer_endpoint, aoi, date_field):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION esri_count_30days STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    url = layer_endpoint.replace('?f=pjson', '') + '/query'

    date = (datetime.utcnow() - timedelta(days=30)).strftime('%Y-%m-%d')

    params = {}
    params['f'] = 'json'
    params['geometryType'] = 'esriGeometryPolygon'
    params['where'] = "{} >= date '{}'".format(date_field, date)
    params['spatialRel'] = 'esriSpatialRelIntersects'
    params['returnGeometry'] = False
    params['returnCountOnly'] = True

    if isinstance(aoi, str):
        aoi = json.loads(aoi)
    count = 0
    for f in aoi['features']:
        params['geometry'] = str({'rings': ogr2rings(f),
                                  'spatialReference': {'wkid': 4326}})
        req = requests.post(url, data=params)
        req.raise_for_status()
        count += req.json()['count']

    # logging.info('FUNCTION esri_count_30days STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return count


def esri_last_instance(layer_endpoint, aoi, field):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION esri_last_instance STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    url = layer_endpoint.replace('?f=pjson', '') + '/query'

    params = {}
    params['f'] = 'json'
    params['geometryType'] = 'esriGeometryPolygon'
    params['where'] = '1=1'
    params['spatialRel'] = 'esriSpatialRelIntersects'
    params['returnGeometry'] = False
    params['outFields'] = field
    params['orderByFields'] = field
    params['returnDistinctValues'] = True

    if isinstance(aoi, str):
        aoi = json.loads(aoi)
    last_instance = None
    for f in aoi['features']:
        params['geometry'] = str({'rings': ogr2rings(f),
                                  'spatialReference': {'wkid': 4326}})
        req = requests.post(url, data=params)
        req.raise_for_status()
        try:
            instances = [item['attributes'][field] for item in
                         req.json()['features']]
            if instances:
                if not last_instance:
                    last_instance = instances[-1]
                elif instances[-1] > last_instance:
                    last_instance = instances[-1]

        except Exception as e:
            raise ValueError((str(e), url, params, req.text))

    # logging.info('FUNCTION esri_last_instance STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return last_instance


# @lru_cache(5)
def cartodb2ogr(service_endpoint, aoi, out_fields, where='', _=''):
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION cartodb2ogr STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    endpoint_template = 'https://{}.carto.com/tables/{}/'
    username, table = search(endpoint_template, service_endpoint + '/')
    url = 'https://{username}.carto.com/api/v2/sql'.format(username=username)

    if isinstance(aoi, str):
        aoi = json.loads(aoi)

    params = {}
    fields = ['ST_AsGeoJSON(the_geom) as geometry']
    out_fields = out_fields.split(',')
    for field in out_fields:
        if field:
            fields.append('{field} as {field}'.format(field=field))

    temp = "ST_Intersects(ST_Buffer(ST_GeomFromText('{}',4326),0),the_geom)"
    features = []
    objectids = []
    for f in aoi['features']:
        where_clause = temp.format(wkt.dumps({'type': 'Polygon',
                                              'coordinates': bbox(f)}))
        if where and not where == '1=1':
            where_clause += 'AND {}'.format(where)

        q = 'SELECT {fields} FROM {table} WHERE {where}'
        params = {'q': q.format(fields=','.join(fields), table=table,
                  where=where_clause)}

        try:
            req = requests.get(url, params=params)
            req.raise_for_status()
        except Exception as e:
            raise ValueError((e, url, bbox(f)))

        response = json.loads(req.text)['rows']
        features += [{
            'type': 'Feature',
            'geometry': json.loads(h['geometry']),
            'properties': {field: h[field] for field in out_fields if field}
        } for h in response]

    featureset = json2ogr({
        'type': 'FeatureCollection',
        'features': features
    })

    # logging.info('FUNCTION cartodb2ogr STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return featureset


def get_split_boxes(f):
    '''
    Check if number of vertices or width or height of bounding box exceed
    thresholds. If they do, returns two revised bounding boxes (Left/Upper
    and Right/Bottom) for intersecting with the geometry
    '''
    x1, y1, x2, y2 = bounds(f)
    if (x2 - x1 > FEATURE_GEOG_THRESHOLD or y2 - y1 > FEATURE_GEOG_THRESHOLD):
        if x2 - x1 > y2 - y1:
            x_split = x1 + (x2 - x1) / 2
            return [box(x1, y1, x_split, y2), box(x_split, y1, x2, y2)]
        else:
            y_split = y1 + (y2 - y1) / 2
            return [box(x1, y1, x2, y_split), box(x1, y_split, x2, y2)]

    return None


def split_multipolygon(f):
    '''
    Split multipolygon into coterminous polygons
    '''
    new_features = [{'type': 'Feature',
                     'properties': f['properties'],
                     'geometry': poly} for poly in f['geometry']]
    return new_features


def split_polygon(f):
    '''
    Split complex geometry in half until they are below vertex and bounding
    box size constraints
    '''
    bbs = get_split_boxes(f)
    new_features = []
    if bbs:
        for bb in bbs:
            geom = f['geometry']
            if not geom.is_valid:
                geom = geom.buffer(0)
            split_feat = {'type': 'Feature',
                        'properties': f['properties'],
                        'geometry': geom.intersection(bb)}
            if split_feat['geometry'].type == 'MultiPolygon':
                poly_feats = split_multipolygon(split_feat)
                for h in poly_feats:
                    new_features.extend(split_polygon(h))
            else:
                new_features.extend(split_polygon(split_feat))
    else:
        new_features.append(f)

    return new_features


def split(featureset):
    '''
    First split all multipolygons into coterminous polygons. Then check each
    against vertex and bounding box size constraints, and split into multiple
    polygons using recursive halving if necessary
    '''
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION split STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    new_features = []
    split_id = 0
    for f in featureset['features']:
        f['properties']['split_id'] = split_id
        split_id += 1
        if f['geometry'].type == 'MultiPolygon':
            poly_feats = split_multipolygon(f)
            for h in poly_feats:
                new_features.extend(split_polygon(h))
        elif f['geometry'].type == 'Polygon':
            new_features.extend(split_polygon(f))

    new_featureset = dict(type=featureset['type'],
                          features=new_features)
    if 'crs' in featureset.keys():
        new_featureset['crs'] = featureset['crs']

    # logging.info('FUNCTION split STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return new_featureset


def get_centroid_bbox(features):
    x, y = zip([f['properties']['centroid'] for f in features])
    return min(x), min(y), max(x), max(y)


def split_features(featureset):
    '''
    Sort features into groups where all features have centroids clustered
    within the geographic size limit and the total number of coordinates is
    beneath the coordinate length limit
    '''
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION split_features STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    feature_groups = [featureset['features']]
    

    # logging.info('FUNCTION split_features STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return


def condense_properties(properties):
    '''
    Combine common properties with duplicate values from all features
    being dissolved
    '''
    return {key: val for key, val in properties[0].items()
            if all(key in p.keys() and val == p[key] for p in properties)}


def dissolve(featureset, fields=None):
    '''
    Dissolve a set of geometries on a field, or dissolve fully to a single
    feature if no field is provided
    '''
    global FUNCTION_COUNT
    FUNCTION_COUNT += 1
    # logging.info('FUNCTION dissolve STEP {} START'.format(FUNCTION_COUNT))
    t0 = time()

    if fields:
        def sort_func(k):
            return ','.join([str(k['properties'][field])
                             for field in fields.split(',')])
    else:
        sort_func = None

    new_features = []
    dissolve_id = 0
    try:
        assert isinstance(featureset, dict)
        assert 'features' in featureset.keys()
        assert isinstance(featureset['features'], list)
    except Exception as e:
        raise ValueError((str(e),featureset))
    if len(featureset['features']) > 0:
        if sort_func:
            features = sorted(featureset['features'], key=sort_func)
            for key, group in itertools.groupby(features, key=sort_func):
                properties, geoms = zip(*[(f['properties'],
                                          f['geometry']) for f in group])
                if geoms and not any(geom is None for geom in geoms):
                    try:
                        new_geom = unary_union(geoms)
                    except Exception as e:
                        new_geom = unary_union([geom if geom.is_valid
                                                else geom.buffer(0)
                                                for geom in geoms])
                    new_properties = condense_properties(properties)
                    new_properties['dissolve_id'] = dissolve_id
                    dissolve_id += 1
                    new_features.append(dict(type='Feature',
                                             geometry=new_geom,
                                             properties=new_properties))

        else:
            geoms = [f['geometry'] if f['geometry'].is_valid else
                     f['geometry'].buffer(0) for f in featureset['features']]
            new_properties = condense_properties([f['properties'] for f in
                                                  featureset['features']])
            new_properties['dissolve_id'] = dissolve_id
            dissolve_id += 1
            new_features.append(dict(type='Feature',
                                         geometry=unary_union(geoms),
                                         properties=new_properties))

    new_featureset = dict(type=featureset['type'],
                          features=new_features)
    if 'crs' in featureset.keys():
        new_featureset['crs'] = featureset['crs']

    # logging.info('FUNCTION dissolve STEP {} DONE - {} SECONDS'.format(FUNCTION_COUNT, time()-t0))
    return new_featureset


def index_featureset(featureset):
    '''
    '''
    index = rtree.index.Index()
    for i, f in enumerate(featureset['features']):