from database.models import Site
from landscape.models import LandCover
from geo.models import LandCoverType
from geo.models import Boundary
from django.db import connection
from django.contrib.gis.geos import GEOSGeometry

from os import path
import numpy
import pandas
import geopandas
import fiona
import pyproj
from matplotlib import pyplot
%matplotlib inline

"""
Explanation: Land cover representation
This notebook compares the land cover within the entire study area to the land cover within a certain radius around sample locations.
Import statements
End of explanation
"""


austria_mgd = "+init=epsg:31254"

"""
Explanation: define the project coordinate system
End of explanation
"""


output_filepath = "/Users/Jake/OneDrive/Documents/alpine soundscapes/data/landcover"

"""
Explanation: define the output filepath
End of explanation
"""


schema = {'geometry': "MultiPolygon", 
          'properties': {'id': 'int', 
                         'cover_type': 'int'}}

"""
Explanation: define schema for GeoJSON output
End of explanation
"""


# returns a geopandas dataframe from a given PostGIS query
def get_geodataframe(queryset, modification=None, index_col=None, crs=austria_mgd):
    query = queryset.query.sql_with_params()
    if modification:
        query = (modification, query[1])
    return geopandas.read_postgis(query[0], connection, 
                                   geom_col='geometry', 
                                   params=query[1], 
                                   index_col=index_col,
                                   crs=crs)

# returns a geopandas dataframe containing the land cover intersected with a given geometry
def get_landcover(geometry, make_valid=False, index_col=None, crs=austria_mgd):
    
    if make_valid:
        intersection = """ST_Intersection(ST_MakeValid("landscape_landcover"."geometry"), %s)"""
    else:
        intersection = """ST_Intersection("landscape_landcover"."geometry", %s)"""
    
    query = """
        SELECT "landscape_landcover"."id", "landscape_landcover"."cover_type", 
        {0} AS "geometry" FROM "landscape_landcover" 
        WHERE ST_Intersects("landscape_landcover"."geometry", %s)
    """.replace('\n', '').format(intersection)

    return geopandas.read_postgis(query, connection, 
                                   geom_col='geometry', 
                                   params=["SRID={0};{1}".format(geometry.srid, geometry.wkt), 
                                           "SRID={0};{1}".format(geometry.srid, geometry.wkt)], 
                                   index_col=index_col,
                                   crs=crs)

# returns a numpy array with the percentage of each land cover type in respect to the total geometry area
def get_percentages(source):
    y = range(1, 16)
    x = numpy.empty(15)
    total_area = source.area.sum()
    for i, y_i in enumerate(y):
        x[i] = (source[source.cover_type == y_i].area.sum() / total_area) * 100
    return x

"""
Explanation: Function declarations
End of explanation
"""


radii = [50, 100, 200, 500]

"""
Explanation: Get location land cover
define a list of radii defining the surrounding land cover
End of explanation
"""


location_area_list = []
for radius in radii:
    location_areas = get_geodataframe(Site.objects.filter(id__lte=30)).buffer(radius)
    location_area = location_areas.iloc[0]
    for area in location_areas:
        location_area = location_area.union(area)
    # convert to geos multipolygon with srid
    location_area = GEOSGeometry(location_area.wkt)
    location_area.srid = 31254
    location_area_list.append(location_area)

"""
Explanation: compute a geometry (unioned buffer around locations) for each radius
End of explanation
"""


location_landcover_list = []
for location_area, radius in zip(location_area_list, radii):
    location_landcover = get_landcover(location_area)
    location_landcover_list.append(location_landcover)
    
    # export result to a GeoJSON file
    #location_landcover.to_file(filename=path.join(output_filepath, "location_landcover_{0}m.geojson".format(radius)), 
    #                           driver='GeoJSON', schema=schema)

"""
Explanation: query database
End of explanation
"""


boundary = Boundary.objects.get(name="study area")

"""
Explanation: Get study area land cover
get the study area boundary geometry from the database
End of explanation
"""


study_landcover = get_landcover(boundary.geometry, make_valid=True)

# export result to a GeoJSON file
#study_landcover.to_file(filename=path.join(output_filepath, "study_area_landcover.geojson".format(radius)), 
#                        driver='GeoJSON', schema=schema)

"""
Explanation: query database
End of explanation
"""


area_percentages = pandas.DataFrame({'name': [r['name'].lower() for r in LandCoverType.objects.filter(id__lte=15).values()],
                                     'study_area': get_percentages(study_landcover)}, index=range(1, 16))

for i, radius in enumerate(radii):
    area_percentages["{0}m".format(radius)] = get_percentages(location_landcover_list[i])

area_percentages

area_percentages.to_csv(path.join(output_filepath, "area_percentages.csv"))

"""
Explanation: Calculate area percentages
End of explanation
"""


