#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 14 13:21:12 2021

Export the stationlist to the geopandas world shape file

@author: geopeter
"""
import numpy as np
import geopandas as gpd


def buildAverageTimeseries(stationList, fromYear, toYear, name):
    meanAverage = []
    for stationID, station in stationList.iterrows():
        temps = []
        for i in range(fromYear,toYear):
            if not np.isnan(station[i]):
                temps.append(station[i])
        
        if len(temps) > 5:
            meanAverage.append(np.mean(temps))
        else:
            meanAverage.append(np.NaN)
    
    stationList[name] = np.round(meanAverage,1)


def cleanAverageTimeseries(stationList):
    # determine gauges that includes both timeseries. If not delete them.
    
    for stationID, station in stationList.iterrows():
        if np.isnan(station['m1961T1990']) or np.isnan(station['m1991T2018']):
            #station['m1961T1990'] = None
            #station['m2010T2018'] = None
            stationList.at[stationID, "m1961T1990"] = None
            stationList.at[stationID, "m1991T2018"] = None


def convertStationListToGPD(stationList):
    
    print("convert stationlist to GeoPandas")
    
    del stationList["file"]
    
    stationList = stationList[stationList.lat != ""]
    stationList = stationList[stationList.lon != ""]
    
    stationList["lat"] = stationList["lat"].astype(str).astype(float)
    stationList["lon"] = stationList["lon"].astype(str).astype(float)
    
    # add some average timeseries
    buildAverageTimeseries(stationList, 1961,1990,"m1961T1990")
    buildAverageTimeseries(stationList, 1991,2018,"m1991T2018")
    cleanAverageTimeseries(stationList)
    
    # convert stationlist to geodataframe
    stationGPD = gpd.GeoDataFrame(stationList, geometry=gpd.points_from_xy(stationList.lat, stationList.lon)).reset_index()
    # del stationGPD["lat"]
    # del stationGPD["lon"]
    stationGPD = stationGPD.set_crs('epsg:4326')
    
    stationGPD.columns = stationGPD.columns.astype(str)
    stationGPD = stationGPD.sort_index(axis=1, ascending=False)
    #stationGPD.to_file("stationList.shp", "ESRI Shapefile")
    
    return stationGPD
    
    
def exportStationListToWorldShape(stationGPD):
    
    print("export stationlist to world shape file")
    
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    countryMeanGPD = gpd.sjoin(world, stationGPD, how="inner", op='intersects')
    countryMeanGPD = countryMeanGPD.groupby("name", axis=0).mean().reset_index()
    
    countryMeanGPD["anomalie"] = countryMeanGPD["m1991T2018"] - countryMeanGPD["m1961T1990"]
    
    del countryMeanGPD["pop_est"]
    del countryMeanGPD["gdp_md_est"]
    
    for i in range(1873,1950):
        try:
            del countryMeanGPD[str(i)]
        except:
            # do nothing
            continue
    
    worldGauge = world.set_index("name").join(countryMeanGPD.set_index("name"))
    worldGauge.to_file("./output/countryAnnualTemperature.shp", "ESRI Shapefile")
    
def export(stationList):
    stationGPD = convertStationListToGPD(stationList)
    exportStationListToWorldShape(stationGPD)