<jupyter_start><jupyter_text>Example of working with GHCN data from Romania[Global historcal climate network](http://www.ncdc.noaa.gov/ghcnm/) weather station data from Romania.Requires Pandas version > 0.19.0<jupyter_code># Load some libraries
import os
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from tqdm import *
from colorama import Fore
matplotlib.style.use('ggplot')
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Step 1: data preperationFrom our Station Data files we need to create:* one single data structure* Date indexed data (with only one index for all the datasets)* Station names as column identifiersWe will use the [Pandas library](http://pandas.pydata.org/pandas-docs/stable/index.html), as it is perfectly suited to our task. We have read it in above with the alias **pd**.Before we try and read all the data, let's test a procedure with one single station file.<jupyter_code># Read a station data file with Pandas
test_data = pd.read_csv("Data/station_data/BUM00015502_VIDIN_BU_.csv")
test_data.head()<jupyter_output><empty_output><jupyter_text>Looks good, but the Dates should be an index, not a column, and they should also be a date object, not a simple integer (we get much more functionality that way).<jupyter_code># Make a list of datetime values out of the integer dates using a list comprehension technique
dates = []
for date in test_data['DATE']:
    dates.append(pd.datetime.strptime(str(date),"%Y%m%d"))
# The above could have been done more effectivley using list comprehension

# Next set the new list as an index, and remove the old column from the dataset
test_data.index = dates
test_data = test_data.drop(['DATE','PRCP'], axis=1)

test_data.head()<jupyter_output><empty_output><jupyter_text>Great! We can plot a simple preview of the data to make sure it looks good.<jupyter_code>test_data.plot()<jupyter_output><empty_output><jupyter_text>Now we have a time-aware time series, we can do some really cool and intelligent sub-scripting of the data:<jupyter_code># Call one day

test_data.TAVG['1998-01-01']
# Call a few days

test_data.TAVG['1998-01-06':'1998-02-07']
# Call a whole month

test_data.TAVG['2000-01']
# Call a large, date range just specifying years

test_data.TAVG['1994':'1996']<jupyter_output><empty_output><jupyter_text>Read mutliple datasetsThe preview plot is messy as our data are not contiguous, but as a quick-check, it seems like everything is more-or-less fine.So, reading a single file is easy, and straightforward. But we want to do some exploratory analysis on multiple station measurements. For this we will need to read all the station data together into a consistent data object.<jupyter_code># Make a small tools (functions) to help with the work

def station_name(fname):
    """Return the station ID from a path/filename.csv string"""
    tmp = fname.split('/')[-1]
    return tmp.split('_')[0]
# If we were on a Mac or Linux system, we could get the file list via a bash command
flist = !ls Data/station_data/*.csv
# But this will break on windows. To make our code cross-platform we use a python
# library to find all the files instead. This is much better than hard-coding the files!

frames = [] # an empty list to hold each data object as it is loaded

mypath = 'Data/station_data/'          # Set path to data
for item in tqdm(os.listdir(mypath)):        # Find all files in that path and loop over them
    if '.csv' in item:                 # If the file is a csv type do something...
        fname = ''.join([mypath,item])
        station = station_name(fname)
        #print('Reading data from station', station)
        tmp = pd.read_csv(fname)
        dates = [pd.datetime.strptime(str(date),"%Y%m%d") for date in tmp['DATE']]
        tmp.index = dates
        tmp = tmp.drop(['DATE','PRCP'], axis=1) # get rid of date and precipitation columns
        tmp.columns = [station]     # Re-name TAVG to be the station name
        frames.append(tmp)
print("{0} GHCN files read".format(len(frames)))
df = pd.concat(frames, axis=1)  # Join all the seperate data together into one object<jupyter_output><empty_output><jupyter_text>Step 2: Cleaning the dataset for analysisNow we have created a dataframe **df** holding all the station data with one coherant time index.This abstraction will do much of the work for us...<jupyter_code># First lets see how long these data run for in time
print("minimum date:", min(df.index).date())
print("maximum date:", max(df.index).date())
# Now let's look at a statistical description of these data
df.describe()<jupyter_output><empty_output><jupyter_text>There is a clear problem with these stats. Most of these data seem to have a missing value of `-9999.0` included.To proceede we should replace with with a missing data type that we can operate with `np.nan`<jupyter_code># we can replace all -999.0 values with np.nan like this
df[df == -9999.0] = np.nan
# Now, the dataframe values seems reasonable, except we can see there are many series which are empty.
# They were just full of missing values for whatever reason.

df.describe()
# It looks like we can simply filter out data that now has a low count (e.g. < 10,000).

limit = 0

for key in df:
    if df[key].count() <= limit:
        print('removing', key,'from df object.')
        df = df.drop([key], axis=1)
# Much better! Finally a clean df object, that we can work from.

df.describe()<jupyter_output><empty_output><jupyter_text>Step 3. Creating a mean time series<jupyter_code>for key in df:
    plt.plot(df[key], lw=0.3, alpha=0.75)
plt.title('Individual stations')
plt.show()<jupyter_output><empty_output><jupyter_text>Aggregated, these data gives a clearer picture<jupyter_code>df.mean(axis=1).plot(title='Romanian TAVG', lw=0.5, alpha=0.8)<jupyter_output><empty_output><jupyter_text>Looking at the data it is easy to see the seasonal signal is dominant. One method of removing the seasonal signal and work with anomalies is to subtract a day-of-the-year (DOY) mean.<jupyter_code># make the average data a new series, and strip out nan values for working ease
df_mean = pd.DataFrame(df.mean(axis=1), columns=['mean_temp'])  # make a new df object
#df_mean = df_mean[df_mean.notnull().values]                     # remove missing values<jupyter_output><empty_output><jupyter_text>I can look at every day in the time series, find its DOY value (an integer from 1 - 366), make a corresponding list, and use it to subscript the mean dataframe, and find the mean on each day of the year. I place these values in a dictionary with integer doy as a key, for ease of use later when deseasonalising.<jupyter_code>doy_index = np.array([ date.dayofyear for date in df_mean.index])  # make an array of DOY values to use as a lookup
d = {}
for doy in range(1,367):
    d[doy] = df_mean[doy_index == doy].mean().values[0]
# Veryify it looks okay...
tmp = []
for key in d:
    tmp.append(d[key])
plt.plot(range(1,367), tmp)
plt.title("Seasonal Temperature signal in Romania")
plt.xlabel('DOY')
plt.xlim([1,367])
plt.ylabel("Temp. °C")
plt.show()<jupyter_output><empty_output><jupyter_text>Step 4: Remove the seasonal signal from these data<jupyter_code>dszn = []
for day in df_mean.index:
    dszn.append(df_mean['mean_temp'][day] - d[day.dayofyear])
dszn = np.array(dszn)

df_mean['anom'] = dszn
# All of the above could have just been done in one line...
# df_mean['anom'] = np.array([df_mean['mean_temp'][day] - d[day.dayofyear] for day in df_mean.index])
plt.plot(df_mean.index, df_mean.anom, lw=0.5, alpha = 0.75)
plt.title(r"Romanina $\delta$ Temp. °C")
plt.show()<jupyter_output><empty_output><jupyter_text>Step 5: Using the data for something useful!Based on historical Romanian average temperature anomalies, how does a given value rank?Requires an average temperature and a date as input.<jupyter_code>def anom_of_given_date(date, temperature):
    """Given a specific date and temperature, return the delta temp value"""
    qtmp =pd.DataFrame([temperature], index=[pd.datetime.strptime(str(date),"%Y%m%d")])
    qanom = temperature - d[qtmp.index[0].dayofyear]
    print("Day of year for {0} is {1}".format(date, qtmp.index[0].dayofyear))
    print("𝛿 temp. = {0:3.2f}°C".format(qanom))
    return qanom
anom = anom_of_given_date(date=19830901, temperature=10)
a, bins, c = plt.hist(df_mean.anom, bins=80, normed=True)
plt.title(r"Romanian $\delta$ temperature")
plt.xlabel(r"$\delta$ Temp. °C")
plt.ylim(0,0.13)
plt.vlines(anom, 0, 0.16, label="day to check")

y = matplotlib.mlab.normpdf(bins, df_mean.anom.mean(), df_mean.anom.std())
plt.plot(bins, y, label='Gaussian pdf')

plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Let's see that in a cumulative density plot as well<jupyter_code>a,b,c = plt.hist(df_mean.anom, bins=100, cumulative=True, normed=True)
plt.vlines(anom, 0, 1.1)
plt.ylim(0,1)
plt.title(r"Cumulative $\delta$ distribution and value to test")
plt.show()<jupyter_output><empty_output><jupyter_text>Rank the anomaly values, and count how many times our value to check is larger than the ranked values. Divide this by the total number of values to get the percentile position.<jupyter_code>ranked = np.array(sorted(df_mean.anom.values))
print("𝛿 of {0:3.2f}°C has a percentile rank of {1:3.2f}".format(anom, len(ranked[anom > ranked]) / len(ranked)))<jupyter_output><empty_output><jupyter_text>Finally, combine the above into a function that can be re-run easily<jupyter_code>def check_anomaly(date, temperature):
    anom = anom_of_given_date(date=date, temperature=temperature)
    a,b,c = plt.hist(df_mean.anom, bins=100, cumulative=True, normed=True)
    plt.vlines(anom, 0, 1.1)
    plt.ylim(0,1)
    plt.title(r"Cumulative $\delta$ distribution")
    plt.ylabel('Density')
    plt.xlabel(r"$\delta$ T. [°C]")
    plt.show()
    ranked = np.array(sorted(df_mean.anom.values))
    print("𝛿 of {0:3.2f}°C has a ".format(anom),end="")
    print(Fore.RED + "percentile rank of {0:3.2f}".format(len(ranked[anom > ranked]) / len(ranked)))
check_anomaly(date=20161020, temperature=15)<jupyter_output><empty_output><jupyter_text>Long-term trends Time resamplingSince we created a data structure that has time as an index, we have the benfit of easily re-sampling (binning correctly through time).E.g. to calculate the annual mean anomaly and trend we can do the following:<jupyter_code># simple way to re-bin data in Pandas ("A") means annually

df_mean.anom.resample("A").mean().plot()<jupyter_output><empty_output><jupyter_text>Let's look at these data with a bit more sophistication (trend, and significance)linear regression taken from [Scipy library](http://www.scipy-lectures.org). Requires two terms: x = time index, y = temperature values<jupyter_code>from scipy import stats
anom_annual = df_mean.anom.resample("A").mean()<jupyter_output><empty_output><jupyter_text>Calculate uncertainy of mean based on [Mean Standard error](https://en.wikipedia.org/wiki/Standard_error).<jupyter_code>sem_annual = df_mean.anom.resample("A").std() / np.sqrt(df_mean.anom.resample("A").count() -1)
slope, intercept, rval, pval, stderr = stats.linregress(x=range(len(anom_annual.values)), y=anom_annual.values)
def my_slope(x, slope, intercept, stderr=None, kind=None):
    """
    Return y-value for a given x based on slope and intercept.
    Changes to return upper or lower bounds, based on error, if keyword is given.
    kind: 'pos' = return upper error bound
          'neg' = return lower error bound
           None   = return mean regression line
    """
    if kind:
        assert stderr, "You forgot to include a value for stderr"
        if kind is 'pos':
            return ((slope * x) + intercept) + stderr
        elif kind is 'neg':
            return ((slope * x) + intercept) - stderr
    else:
        return (slope * x) + intercept
# Use the slope, intercept and error to calculate a linear fit to the data

fit = [ my_slope(x=x, slope=slope, intercept=intercept) for x in range(len(anom_annual)) ]
fit_pos = [ my_slope(x=x, slope=slope, intercept=intercept, stderr=stderr, kind='pos') for x in range(len(anom_annual)) ]
fit_neg = [ my_slope(x=x, slope=slope, intercept=intercept, stderr=stderr, kind='neg') for x in range(len(anom_annual)) ]
# Plot everything with errors, and print a summary of the stats...

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

ax.errorbar(anom_annual.index, anom_annual.values, yerr=sem_annual.values)

ax.plot(anom_annual.index, fit, color='red', lw=0.75)
ax.plot(anom_annual.index, fit_pos, 'r--', lw=0.75)
ax.plot(anom_annual.index, fit_neg, 'r--', lw=0.75)

#fig.title("Romania GHCN Annual temp. anomaly")
#ax.ylabel("$\delta$ Temp. °C")
#ax.xlabel("Year")

my_text = r"$y$ = {0:3.2f} $\times x +$ {1:3.2f}$\pm$ {2:3.2f}".format(slope, intercept, stderr)
my_text += '\n'+r'$r^2={0:3.2f}, p=${1:3.5f}'.format(rval*rval, pval)

fig.text(0.15, 0.8, s=my_text)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12036535335552287, 'fraction_numerical': 0.016564749593621798, 'mean_word_length': 4.724412937527692, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 84, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '398173', 'license': 'MIT', 'n_tokens_mistral': 4532, 'n_tokens_neox': 4327, 'n_words': 1414, 'path': 'test_stations.ipynb', 'repo_name': 'benlaken/Romania_2016', 'size': '12919'}