<jupyter_start><jupyter_text>Resample to 24h accumulation<jupyter_code>charlotte_24h_rainfall = pd.DataFrame()
charlotte_24h_rainfall['mean_rain_1'] = charlotte_rainfall.Rainfall_1.resample('D').mean()
charlotte_24h_rainfall['accum_rain_1'] = charlotte_rainfall.Rainfall_1.resample('D').sum()
charlotte_24h_rainfall['mean_rain_2'] = charlotte_rainfall.Rainfall_2.resample('D').mean()
charlotte_24h_rainfall['accum_rain_2'] = charlotte_rainfall.Rainfall_2.resample('D').sum()
mask = (charlotte_24h_rainfall.accum_rain_1 == 0) | (charlotte_24h_rainfall.accum_rain_2 == 0)
charlotte_24h_rainfall = charlotte_24h_rainfall.loc[~mask]
mask = (np.isnan(charlotte_24h_rainfall.accum_rain_1)) | (np.isnan(charlotte_24h_rainfall.accum_rain_2))
charlotte_24h_rainfall = charlotte_24h_rainfall.loc[~mask]
fig = plt.figure()
plt.plot(charlotte_24h_rainfall.index, charlotte_24h_rainfall['accum_rain_1'])
plt.ylabel('mm/24h')
fig.autofmt_xdate()
fig = plt.figure()
plt.plot(charlotte_24h_rainfall.index, charlotte_24h_rainfall['accum_rain_2'])
plt.ylabel('mm/24h')
fig.autofmt_xdate()<jupyter_output><empty_output><jupyter_text>Try plotting Semivariogram<jupyter_code>#charlotte_24h_rainfall['minutes'] = np.arange(0, 24*60*len(charlotte_24h_rainfall.accum_rain_1), 24*60)
charlotte_24h_rainfall['seconds'] = charlotte_24h_rainfall.index.astype(np.int64) // 10 ** 9 - ( charlotte_24h_rainfall.index.astype(np.int64)[0] // 10 ** 9)<jupyter_output><empty_output><jupyter_text>Set seconds from zero to length<jupyter_code>charlotte_24h_rainfall['seconds'] = np.arange(0,len(charlotte_24h_rainfall.accum_rain_1), 1)
charlotte_24h_rainfall['zeros'] = np.zeros((len(charlotte_24h_rainfall.accum_rain_1), 1), dtype=np.int8)
charlotte_24h_rainfall.head()<jupyter_output><empty_output><jupyter_text>New try<jupyter_code>from rpy2.robjects.packages import importr
from rpy2.robjects import r
sp = importr('sp')
gstat = importr('gstat')
intamap = importr('intamap')
r('jet.colors  <- c("#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000")')
r('col.palette <- colorRampPalette(jet.colors)')
rain1 = charlotte_24h_rainfall[['accum_rain_1', 'seconds', 'zeros']]
rain2 = charlotte_24h_rainfall[['accum_rain_2', 'seconds', 'zeros']]
from rpy2.robjects import pandas2ri
pandas2ri.activate()
r_df = pandas2ri.py2ri(rain1)
r.assign('mydata', r_df)
r_df2 = pandas2ri.py2ri(rain2)
r.assign('mydata2', r_df2)
r('''
mydata <- data.frame(mydata)
coordinates(mydata) <- ~seconds+zeros

mydata2 <- data.frame(mydata2)
coordinates(mydata2) <- ~seconds+zeros
''')
p_myiso = r('myiso <- variogram(accum_rain_1~1,mydata,width=5,cutoff=100)')
p_myiso2 = r('myiso2 <- variogram(accum_rain_2~1,mydata2,width=5,cutoff=100)')
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(p_myiso['dist'], p_myiso['gamma'], '-o')
plt.plot(p_myiso2['dist'], p_myiso2['gamma'], '-o')
p_myiso
p_myiso2
r('''
variogram_map <- variogram(accum_rain_1~1,mydata,width=1,cutoff=50,map=TRUE)
png("variogram_map_temporal.png",height=600,width=600)
print(plot(variogram_map))
dev.off()
''')
r('''
initial_vario_sph <- vgm(psill=100,model="Lin",range=8,nugget=350)
fitted_vario_sph <- fit.variogram(myiso,initial_vario_sph)
range  <- fitted_vario_sph$range[2]				# fitted range
nugget <- fitted_vario_sph$psill[1]				# fitted nugget
sill   <- sum(fitted_vario_sph$psill)			# fitted sill (total sill)
SSErr_sph <- attributes(fitted_vario_sph)$SSErr	# sum of squared errors (goodness of fit)

png("fitted_isotropic_variogram_sph_temporal.png",height=600,width=900)
print(plot(myiso,fitted_vario_sph))
dev.off()
''')
r('range')
r('nugget')
r('sill')
r('SSErr_sph')<jupyter_output><empty_output><jupyter_text>Rainfall values on 5th august 2011<jupyter_code>event_rain = charlotte_rainfall.loc['2011-08-05 00:00:00':'2011-08-05 23:59:59']
event_rain_1 = event_rain.Rainfall_1
event_rain_2 = event_rain.Rainfall_2
import datetime as dt
event_rain_2.index = event_rain_2.index + dt.timedelta(hours=72)
event_rain_1 = pd.Series(event_rain_1.Rainfall_1)
my_event = event_rain_1.append(event_rain_2)
my_event = pd.DataFrame(my_event)
my_event.columns = ['R']

my_event['seconds'] = my_event.index.astype(np.int64) // 10 ** 9 - ( my_event.index.astype(np.int64)[0] // 10 ** 9)
my_event['zeros'] = np.zeros((len(my_event.R), 1), dtype=np.int8)
my_event.tail()
r_df = pandas2ri.py2ri(my_event)
r.assign('mydata_event', r_df)
r('''
mydata_event <- data.frame(mydata_event)
coordinates(mydata_event) <- ~seconds+zeros
''')
p_myiso = r('myiso_event <- variogram(R~1,mydata_event,width=1800,cutoff=250000)')
plt.plot(p_myiso['dist'], p_myiso['gamma'], '-o')
from pandas.plotting import autocorrelation_plot
autocorrelation_plot(my_event.R)
autocorrelation_plot(rain1.accum_rain_1)
autocorrelation_plot(rain2.accum_rain_2)
from scipy.spatial.distance import pdist, squareform
data = np.array(charlotte_24h_rainfall[['seconds', 'zeros', 'accum_rain_1']])
lags = np.arange(86400, 864000, 86400)
sill = np.var(data[:,2])
print(sill)
print(lags)
distances = squareform(pdist(data[:,:2]))
print(distances)
tail = data[0:-1,2]
head = data[1::,2]
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum
print(cur_var)
tail = data[0:-2,2]
head = data[2::,2]
print(len(head))
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum

print(cur_var)
tail = data[0:-3,2]
head = data[3::,2]
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum
print(cur_var)
tail = data[0:-4,2]
head = data[4::,2]
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum
print(cur_var)
tail = data[0:-5,2]
head = data[5::,2]
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum
print(cur_var)
tail = data[0:-6,2]
head = data[6::,2]
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum
print(cur_var)
tail = data[0:-7,2]
head = data[7::,2]
cur_diff = (tail-head)**2
cur_sum = sum(cur_diff)
cur_var = 1/(2*len(tail))*cur_sum
print(cur_var)
data = np.array(charlotte_24h_rainfall[['seconds', 'zeros', 'accum_rain_1']])
npoints, cols = data.shape
npoints
cols
from scipy.spatial.distance import pdist, squareform
square_data = squareform(pdist(data[:,:2]))
lags = np.arange(1*7*24*60*60, 6*30*24*60*60, 1*7*24*60*60 )
sill = np.var(data[:,2])
sill
lags.shape
def lagindices( pwdist, lag, tol ):
    '''
    Input:  (pwdist) square NumPy array of pairwise distances
            (lag)    the distance, h, between points
            (tol)    the tolerance we are comfortable with around (lag)
    Output: (ind)    list of tuples; the first element is the row of
                     (data) for one point, the second element is the row
                     of a point (lag)+/-(tol) away from the first point,
                     e.g., (3,5) corresponds fo data[3,:], and data[5,:]
    '''
    # grab the coordinates in a given range: lag +/- tolerance
    i, j = np.where( ( pwdist >= lag - tol )&( pwdist < lag + tol ) )
    # zip the coordinates into a list
    indices = zip( i, j )
    # take out the repeated elements,
    # since p is a *symmetric* distance matrix
    indices = np.array([ i for i in indices if i[1] > i[0] ])
    return indices
def semivariance( data, indices ):
    '''
    Input:  (data)    NumPy array where the fris t two columns
                      are the spatial coordinates, x and y, and
                      the third column is the variable of interest
            (indices) indices of paired data points in (data)
    Output:  (z)      semivariance value at lag (h) +/- (tol)
    '''
    # take the squared difference between
    # the values of the variable of interest
    z = [ ( data[i,2] - data[j,2] )**2.0 for i,j in indices ]
    # the semivariance is half the mean squared difference
    return np.mean( z ) / 2.0
index = [ lagindices( square_data, lag, 0.1 ) for lag in lags ]
lags
v = np.array([ semivariance( data, indices ) for indices in index ])
v
v.size
vario = np.array( zip( lags, v ) ).T
h, sv = np.array( zip( lags, v ) ).T
h, sv = lags, v
plt.plot( h, sv, 'ko-' )
def SVh( P, h, bw ):
    '''
    Experimental semivariogram for a single lag
    '''
    pd = squareform( pdist( P[:,:2] ) )
    N = pd.shape[0]
    Z = list()
    for i in range(N):
        for j in range(i+1,N):
            if( pd[i,j] >= h-bw )and( pd[i,j] <= h+bw ):
                Z.append( ( P[i,2] - P[j,2] )**2.0 )
    return np.sum( Z ) / ( 2.0 * len( Z ) )
 
def SV( P, hs, bw ):
    '''
    Experimental variogram for a collection of lags
    '''
    sv = list()
    for h in hs:
        sv.append( SVh( P, h, bw ) )
    sv = [ [ hs[i], sv[i] ] for i in range( len( hs ) ) if sv[i] > 0 ]
    return np.array( sv ).T
 
def C( P, h, bw ):
    '''
    Calculate the sill
    '''
    c0 = np.var( P[:,2] )
    if h == 0:
        return c0
    return c0 - SVh( P, h, bw )
86400/2.0
lags
lags = np.arange(86400, 86400*15, 86400/2.0)
print(lags)
sv = SV( data, lags, 10000 )
plt.plot( sv[0], sv[1], '.-' )<jupyter_output>/home/kriete/.virtualenvs/cie5703/lib/python3.5/site-packages/ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in double_scalars
  if sys.path[0] == '':
{'fraction_non_alphanumeric': 0.16116061725711803, 'fraction_numerical': 0.04825038035209737, 'mean_word_length': 4.861783439490446, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '255144', 'license': 'MIT', 'n_tokens_mistral': 4118, 'n_tokens_neox': 3716, 'n_words': 719, 'path': 'week_7_variograms.ipynb', 'repo_name': 'kriete/cie5703_notebooks', 'size': '9202'}