stations = pd.read_csv(
    '../../data/ncdc/ish-history.csv',
    skiprows=1, 
    names=['usaf','wban','stname','ctry','fips','state','call','lat' ,'lon' ,'elev'],
    dtype=object)
print (len(stations))
stations[:3]

"""
Explanation: Import stations CSV
My first step is to find interesting stations and one of the requirements is to find some that have a good time range coverage so I will evaluate the files downloaded.
Reading the stations coordinates CSV
End of explanation
"""


stations['lat'] = stations.lat.apply(lambda lat: float(lat)/1000)
stations['lon'] = stations.lon.apply(lambda lon: float(lon)/1000)
stations['elev'] = stations.elev.apply(lambda e: float(e)/10.0)

stations['id']     = stations.apply(lambda row : "{}-{}".format(row.usaf,row.wban),axis=1)
stations.set_index(['id'],inplace=True)
stations.dropna(how='any',subset=['lat','lon'],inplace=True)
stations.head()

stations[(stations.usaf == '082840')]

"""
Explanation: Get correct coordinates dividing the lat/lon by 1000 and the elevation by 10. Then generate a unique index and a column that specficies if the id is from a USAF station or from WBAN. Finally drop any station that doesn't have a location to use.
End of explanation
"""


p = Path('../../data/ncdc')

"""
Explanation: Read CSV of observations
Instantiate the path where the NCDC data has been downloaded
End of explanation
"""


gsodCSV = p.joinpath('gsod.csv')
if not gsodCSV.exists():
    ops = p.joinpath('raw').joinpath('gsod').glob('**/*.op')
    agsod = []

    for op in ops:
        try:
            data = op.name.replace('.op','').split('-')
            data.append(sum(1 for line in op.open( encoding='utf8' ))-1)
            agsod.append(data)
        except UnicodeDecodeError:
          print (op.absolute())
    dfGsod = pd.DataFrame(data = agsod, columns=['usaf','wban','year','obs'])
    dfGsod['id'] = dfGsod.apply(lambda row: "{}-{}".format(row.usaf,row.wban) ,axis=1)
    dfGsod = dfGsod.set_index(['id'])
    dfGsod[['year','obs']].to_csv(str(gsodCSV))



print ('Reading existing stations per year CSV')
dfGsod = pd.read_csv(str(gsodCSV),index_col=['id'])
print ("{:,} station files".format(len(dfGsod)))
print ("{:,} observations".format(dfGsod.obs.sum()))
dfGsod.head()

"""
Explanation: Check if we have already the stations CSV or generate it from the files. To generate the CSV the name file will be inspected but also a full read of the file will be performed in order to count the number of observations on every station.
End of explanation
"""


year_groups = dfGsod[['year']].groupby(level=0)
year_count  = year_groups.count()
year_max    = year_groups.max()
year_min    = year_groups.min()

years = pd.concat([year_count,year_max,year_min],axis=1)
years.columns = ['count','max','min']
years.head()

"""
Explanation: Year statistics per station
Now study this dataframe, grouping by id and see the total years recorded, max and min
End of explanation
"""


obs_groups = dfGsod[['obs']].groupby(level=0)
obs_count = obs_groups.sum()
obs_count.head()

obs_count[(obs_count.index == '082840-99999')]

"""
Explanation: Count observations per station
End of explanation
"""


stations.index.is_unique and years.index.is_unique and obs_count.index.is_unique

scdf = pd.concat([stations[(stations.index=='082840-99999')],years,obs_count],axis=1,join='inner')
scdf.head()

"""
Explanation: Join stations and observations statistics
Now we can check if the indexes of both data frames are unique and then join them to retreive only the stations with observations
End of explanation
"""


scdf.to_csv('stations_vlc.csv')

"""
Explanation: Finally we can study this dataset and filter the appropriate stations for our study on this first iteration
End of explanation
"""


selections.index.is_unique

scdf['koppen'] = 'BSk'

"""
Explanation: Check if the index of the imported dataframe is unique and then join it with our stations dataset. This join will keep only data on both data frames using the parameter join='inner'.
End of explanation
"""


files_to_read = pd.merge(left=dfGsod,right=scdf,left_index=True,right_index=True,)[['year']]
files_to_read['id'] = files_to_read.index
files_to_read['path'] = files_to_read.apply(lambda row: Path('../../data/ncdc/raw/gsod').joinpath("{0}/{1}-{0}.op".format(row.year,row.id)),axis=1)
print ("{} files to read".format(len(files_to_read)))

"""
Explanation: Getting the observations for the selected stations
Filter the observations data set (around 11M records) by the selected stations and generate the file list to read from the observations folder. First we join the stations with the counts, then we move the index to a column and finally we define the path as a concatenation of the year (on a folder) and the file name as the id and the year.
End of explanation
"""


def getId(stn,wban):
    try:
        istn = int(stn)
        iwban = int(wban)
        return "{:0>6}-{:0>5}".format(istn,iwban)
    except ValueError:
        print("{}/{}".format(stn,wban))
        
def getStationByStnWban(stn,wban):
    try:
        koppen = scdf.loc[getId(stn,wban)].koppen
    except KeyError:
        koppen = None
    return koppen        

observationsCSV = p.joinpath('observations_vlc.csv')
if not observationsCSV.exists():
    i = 0
    acc = 0
    for index,row in files_to_read.iterrows():
        path = row['path']
        if path.exists():
            dfObsTemp = pd.read_fwf(
                    str(path),
                    colspecs=[(0,7),(7,13),(14,18),(18,22),(25,30),(31,33),(35,41),
                              (42,44),(46,52),(53,55),(57,63),(64,66),(68,73),(74,76),
                              (78,84),(84,86),(88,93),(95,100),(102,108),(108,109),
                              (110,116),(116,117),(118,123),(123,124),(125,130),(132,138)],
                    skiprows=1,
                    names=['stn','wban','year','monthday','temp','temp_count',
                           'dewp','dewp_count','slp','slp_count','stp','stp_count',
                           'visib','visib_count','wsdp','wsdp_count','mxspd',
                           'gust','max','max_flag','min','min_flag','prcp','prc_flag','sndp','frshtt']
                    )
            dfObsTemp['koppen'] = dfObsTemp.apply(lambda row: getStationByStnWban(row.stn,row.wban),axis=1)
            dfObsTemp.to_csv(str(observationsCSV),mode='a')

            i += 1
            acc += len (dfObsTemp)

            if i % 1000 == 0:
                print("{:>8} obs".format(acc))

print ('Reading observations CSV')
dfObs = pd.read_csv(str(observationsCSV),index_col=0)
dfObs = dfObs[(dfObs.stn != 'stn')]
print ("{:,} observations".format(len(dfObs)))
dfObs.head()

"""
Explanation: Read the files defined previously and store the results on a new big data frame and CSV adding the Köppen classification. First we define a couple of helping functions and then we iterate over the files_to_read data frame to read the corresponding CSV and append it to the final CSV file.
End of explanation
"""


dfObs2 = dfObs.copy()[['stn','wban','year','monthday','temp','max','min','frshtt','koppen']]

"""
Explanation: Performing data management operations on the dataset
Now that we have the raw dataset, we can start doing management operations. First generate a copy dataframe with only the columns we are interested in.
End of explanation
"""


def getDateTimeFromRow(row):
    try:
        iyear = int(row.year)
        imonth = int("{:0>4}".format(row.monthday)[0:2])
        iday = int("{:0>4}".format(row.monthday)[2:4])
        return  datetime(iyear,imonth,iday)
    except ValueError:
        return np.nan

dfObs2['id']   = dfObs2.apply(lambda row: getId(row.stn,row.wban) ,axis=1)
dfObs2['date'] = dfObs2.apply(lambda row : getDateTimeFromRow(row),axis=1)
dfObs2.set_index(['id','date'],inplace=True)

"""
Explanation: Management
Generate an index using the id station and the date
End of explanation
"""


dfObs2['frshtt']  = dfObs2.apply(lambda row: "{:0>6}".format(row.frshtt),axis=1)
dfObs2['fog']     = dfObs2['frshtt'].apply(lambda row: row[0:1]=='1')
dfObs2['rain']    = dfObs2['frshtt'].apply(lambda row: row[1:2]=='1')
dfObs2['snow']    = dfObs2['frshtt'].apply(lambda row: row[2:3]=='1')
dfObs2['hail']    = dfObs2['frshtt'].apply(lambda row: row[3:4]=='1')
dfObs2['thunder'] = dfObs2['frshtt'].apply(lambda row: row[4:5]=='1')
dfObs2['tornado'] = dfObs2['frshtt'].apply(lambda row: row[5:6]=='1')

"""
Explanation: The frshtt column needs to be padded with zeros to get all the flags in the correct place. Then is possible to get the occurrence of different weather conditions
End of explanation
"""


dfObs2['tempC'] = dfObs2['temp'].replace('99.9', np.nan)
dfObs2['maxC']  = dfObs2['max'].replace('99.9', np.nan)
dfObs2['minC']  = dfObs2['min'].replace('99.9', np.nan)

dfObs2['tempC'] = dfObs2['tempC'].convert_objects(convert_numeric=True)
dfObs2['maxC']  = dfObs2['maxC'].convert_objects(convert_numeric=True)
dfObs2['minC']  = dfObs2['minC'].convert_objects(convert_numeric=True) 

def FtoC(f):
    return (f-32)*5/9

dfObs2['tempC'] = dfObs2['tempC'].apply(lambda temp: FtoC(temp))
dfObs2['maxC'] = dfObs2['maxC'].apply(lambda temp: FtoC(temp))
dfObs2['minC'] = dfObs2['minC'].apply(lambda temp: FtoC(temp))

dfObsTemp.to_csv(str(observationsCSV))

"""
Explanation: Recode the temperatures columns, replacing the NaN values and afterwards as numerics
End of explanation
"""


dfObs.koppen.value_counts(normalize=True)*100

dfObs2.tornado.value_counts(normalize=True)*100

dfObs2.thunder.value_counts(normalize=True)*100

"""
Explanation: Frequency tables
Frequency tables for koppen, thunders, tornados that are categorized
End of explanation
"""


dfObs2['temp4']=pd.qcut(dfObs2.tempC, 4, labels=["1=0%tile","2=25%tile","3=50%tile","4=75%tile"])
dfObs2['temp4'].value_counts(normalize=True)*100

"""
Explanation: Categorize the temperatures by quantiles and then make the frequency table to confirm the categorization
End of explanation
"""


group_by_year= dfObs2[['temp4','tempC']].groupby(['temp4'])
group_by_year.max()

"""
Explanation: And to get the cuts, we can group by the categorization column and get the max value for every group
End of explanation
"""


