import numpy as np
import matplotlib.pyplot as plt
import scipy.signal

# IPython "magic" to allow plots to be displayed inline
%matplotlib inline

SIZE = 1000
x = np.arange(0, SIZE)

noise = np.random.uniform(low=-50, high=50, size=SIZE)
sin_wave = 100 * np.sin( (np.pi/180) * x )
trend = np.linspace(-500, 500, SIZE)

plt.plot(x, noise, label='noise')
plt.plot(x, sin_wave, label='sin wave')
plt.plot(x, trend, label = 'trend')
plt.legend()

noisy_sin = noise + sin_wave
noisy_trending_sin = noisy_sin + trend

plt.plot(x, noisy_sin, label='noisy sin wave')
plt.plot(x, noisy_trending_sin, label='noisy, trending, sin wave')
plt.legend()

lin_detrend = scipy.signal.detrend(noisy_trending_sin, type='linear')
con_detrend = scipy.signal.detrend(noisy_trending_sin, type='constant')

plt.plot(x, noisy_trending_sin, label='noisy trending sin')
plt.plot(x, lin_detrend, label='linear detrend')
plt.plot(x, scipy.signal.detrend(noisy_trending_sin, type='constant'), label='constant detrend')

plt.legend()

offset_lin_detrend = lin_detrend + (noisy_trending_sin[0] - lin_detrend[0])
plt.plot(x, noisy_trending_sin, label='noisy trending sin')
plt.plot(x, offset_lin_detrend,  label='offset linear detrend')
plt.legend()

%reset -f

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal

import rasterio
import netCDF4  # Not using this yet, but may want to write netcdf file??

# Create a list of the files we'd like to include 
# in our virtual dataset. For starters, we'll 
# just pick January for all years:
!ls ../../snap-data/tas_mean_C_iem_cccma_cgcm3_1_sresa1b_2001_2100/tas_mean_C_iem_cccma_cgcm3_1_sresa1b_01* > janfiles.txt

# Next, we ask gdalbuildvrt to create a .vrt file
# from all the files in the list we just created.
#   -seperate      each file to a stacked band
#   -b 1           use band number 1 from each input file
%time !gdalbuildvrt -separate -b 1 -input_file_list janfiles.txt jan.vrt

month_temperature_file = rasterio.open("jan.vrt")

# this one takes a while...
%time temperature_data = month_temperature_file.read(masked=None)

# dimensions are (time, y, x) 
# from upper left corner of image
# (not sure which corner of the pixel)

print "Type of data: ", type(temperature_data)
print "Size of data: ", temperature_data.nbytes/1024/1024, "(Mbs)"
print "Shape of data: ", temperature_data.shape
print "Min value: ", np.min(temperature_data)
print "Max value: ", np.max(temperature_data)



print "Using mask?:", temperature_data.mask
print "Fill value:", temperature_data.fill_value
print "First value (should be masked):", temperature_data[0,0,0]


temperature_data = np.ma.masked_outside(temperature_data, -500, 500, copy=False)

print "First value: (should be masked):", temperature_data[0,0,0]
print ""
plt.imshow(temperature_data[0,:,:])
plt.colorbar()

# Next, we can look at the timeseries for a 
# pixel in central alaska that should have data
plt.plot(temperature_data[:, 550, 1050])

# Another 15 seconds or so and a lot of memory. There might
# be a better way to do this.
%time dtrd = scipy.signal.detrend(temperature_data, axis=0)

Y = 550
X = 1050
plt.suptitle("Pixel(y,x): (%s, %s)" % (Y, X))
plt.plot(temperature_data[:,Y,X], label='original')
plt.plot(dtrd[:,Y,X], label='detrended')
plt.plot((dtrd[:,Y,X] + (temperature_data[0,Y,X] - dtrd[0,Y,X])), label='dtrend+offset')
plt.legend()

%time dtrd[:,:,:] += (temperature_data[0,:,:] - dtrd[0,:,:])



def check_pixel(orig, detrended, Y = 550, X = 1050):
    plt.suptitle("Pixel(y,x): (%s, %s)" % (Y, X))
    plt.plot(orig[:,Y,X], label='original')
    plt.plot(detrended[:,Y,X], label='detrended')
    plt.legend()

check_pixel(temperature_data, dtrd, Y=550,X=1050)    

check_pixel(temperature_data, dtrd, Y=1200, X=780)

# Yikes, lots of memory use! 
# See what we are using:
print dtrd.nbytes/1024.0/1024.0/1024.0, "GBs"
print temperature_data.nbytes/1024.0/1024.0/1024.0, "GBs"







%reset -f



import numpy as np
import matplotlib.pyplot as plt
import rasterio
ds_orig = rasterio.open('orig.vrt', 'r')
ds_new = rasterio.open('month-12-new.vrt', 'r')
ds_orig_data = ds_orig.read?
ds_orig_data = ds_orig.read()
ds_new_data = ds_new.read()
plt.imshow(ds_orig_data[0])
plt.ion()
plt.show()
ds_orig_data = np.ma.masked_outside(ds_orig_data, -500, 500, copy=False)
ds_new_data = np.ma.masked_outside(ds_new_data, -500, 500, copy=False)
plt.imshow(ds_orig_data[0])
f0 = plt.figure()
f1 = plt.figure()
plt.plot(ds_new_data[0,550,1050])
plt.plot(ds_new_data[:,550, 1050])
plt.plot(ds_orig_data[:,550,1050], label='orig')
plt.plot(ds_new_data[:,550,1050], label='detrended')
ds_orig_data[0,0,0]
ds_orig_data[0,550,1050]
ds_new_data[0,550,1050]
plt.plot(ds_new_data[:,551,1050], label='detrended')
plt.plot(ds_orig_data[:,551,1050], label='orig')
plt.plot(scipy.signal.detrend(ds_orig_data[:,551,1050], label='fresh detrend')
)
import scipy.signal
plt.plot(scipy.signal.detrend(ds_orig_data[:,551,1050], label='fresh detrend')
)
plt.plot(scipy.signal.detrend(ds_orig_data[:,551,1050], axis=0), label='fresh detrend')
)
plt.plot(scipy.signal.detrend(ds_orig_data[:,551,1050], axis=0), label='fresh detrend')
fdt = scipy.signal.detrend(ds_orig_data[:,551,1050], axis=0
)
fdt[0]
ds_orig_data[0,551,1050]
ds_new_data[0,551,1050]
fdt[0]+=ds_orig_data[0,551,1051]
fdt[0]-=ds_orig_data[0,551,1051]
fdt[0] + ds_orig_data[0,551,1051]
plt.plot(fdt + ds_orig_data[0,551,1050], label='fresh detrend + offset')
plt.legend()
plt.plot(fdt + (ds_orig_data[0,551,1050]-fdt[0]), label='fresh detrend + offset2')
plt.legend()
history