import sys,os
# put local directory into the path
sys.path.insert(0,os.path.abspath('files%spython'%os.sep))

# import module
from gzurl import gzurl

help(gzurl)

import urllib2, io, gzip, tempfile

# codes for url specification on globalbedo.org
years = range(1998,2012)
codes = [95,95,97,97,26,66,54,54,29,25,53,56,56,78]
XX = dict(zip(years,codes))

year = 2009

root = 'http://www.globalbedo.org/GlobAlbedo%d/mosaics/%d/0.5/monthly/'%\
        (XX[year],year)
    
# filename formatting string: use %02d for month eg 01 for 1
month = 1
url = root + '/GlobAlbedo.%d%02d.mosaic.5.nc.gz'%(year,month)
print url

# open file from url
f = urllib2.urlopen(url)

bdata = f.read()

# which looks like this:
bdata[:50]

f = urllib2.urlopen(url)
fileobj = io.BytesIO(f.read())

gzip.GzipFile(fileobj=fileobj)

f = urllib2.urlopen(url)
data=gzip.GzipFile(fileobj=io.BytesIO(f.read())).read()

tmp = tempfile.NamedTemporaryFile(delete=False)
print tmp.name

tmp.write(data)

tmp.unlink(tmp.name)

import sys,os
# put local directory into the path
sys.path.insert(0,os.path.abspath('python'))
# import local module gzurl
from gzurl import gzurl
import gdal

def readGA(root='data/',year=2009,month=1,layer = 'BHR_VIS',filename=None):
    '''
    Method to read a GlobAlbedo file from earlier
    '''
    file_template = 'NETCDF:"%s":%s'

    # allow filename to be overridden from filename=
    filename = filename or root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)

    g = gdal.Open (  file_template % ( filename, layer ) )

    if g is None:
      raise IOError
    data = g.ReadAsArray()
    
    # return a numpy array
    return(np.array(data))


# codes for url specification on globalbedo.org
years = range(1998,2012)
codes = [95,95,97,97,26,66,54,54,29,25,53,56,56,78]
XX = dict(zip(years,codes))

year = 2009

root = 'http://www.globalbedo.org/GlobAlbedo%d/mosaics/%d/0.5/monthly/'%\
        (XX[year],year)
    
print root

# filename formatting string: use %02d for month eg 01 for 1
month = 1
url = root + '/GlobAlbedo.%d%02d.mosaic.5.nc.gz'%(year,month)

# read the gzipped file
f = gzurl(url)
# read the netCDF file from f.filename
nc = readGA(filename=f.filename)
print nc

import sys,os
# put local directory into the path
sys.path.insert(0,os.path.abspath('files%spython'%os.sep))
# import local module gzurl
from gzurl import gzurl
import gdal

# codes for url specification on globalbedo.org
years = range(1998,2012)
codes = [95,95,97,97,26,66,54,54,29,25,53,56,56,78]
XX = dict(zip(years,codes))

year = 2009

root = 'http://www.globalbedo.org/GlobAlbedo%d/mosaics/%d/0.5/monthly/'%\
        (XX[year],year)

for month0 in range(12):  
    # filename formatting string: use %02d for month eg 01 for 1
    base = 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month0+1)
    url = root + base + '.gz'
    # specify a local filename
    # work out how / why this works ...
    local = os.path.join('data{0}'.format(os.sep),base)
    
    # read the gzipped file
    print local
    f = gzurl(url,filename=local)
    # read the netCDF file from f.filename
    # read the gzipped file
    nc = readGA(filename=f.filename)

!rm -f data/GlobAlbedo.200901.mosaic.5.nc.gz 
!wget -O data/GlobAlbedo.200901.mosaic.5.nc.gz \
  http://www.globalbedo.org/GlobAlbedo56/mosaics/2009/0.5/monthly/GlobAlbedo.200901.mosaic.5.nc.gz

!gunzip -f data/GlobAlbedo.200901.mosaic.5.nc.gz
!ls -l data/GlobAlbedo.200901.mosaic.5.nc

import sys,os
sys.path.insert(0,os.path.abspath('files%spython'%os.sep))
from gzurl import gzurl
from netCDF4 import Dataset

years = range(1998,2012)
codes = [95,95,97,97,26,66,54,54,29,25,53,56,56,78]
XX = dict(zip(years,codes))

year = 2009
tile = 'h17v03'

root = 'http://www.globalbedo.org/GlobAlbedo%d/tiles/%d/%s/'%\
        (XX[year],year,tile)
    
# filename formatting string: use %03d for doy eg 001 for 1
doy = 145
url = root + 'GlobAlbedo.%d%03d.%s.nc.gz'%(year,doy,tile)
# see if you can make sense of this complicated formatting
filename = url.split('/')[-1].replace('.gz','')
local_file = 'files{0}data{0}{1}'.format(os.sep,filename)

# try to read local file
try:
    nc = Dataset(local_file,'r')
except:
    f = gzurl(url,filename=local_file)
    nc = Dataset(f.filename,'r')
    f.close()


# now pull some data

vis = np.array(nc.variables['BHR_VIS'])
nir = np.array(nc.variables['BHR_NIR'])
ndvi = (nir - vis)/(nir + vis)

import pylab as plt

# figure size
plt.figure(figsize=(8,8))
# title
plt.title('NDVI: Tile %s %d doy %03d'%(tile,year,doy))
# colour map
cmap = plt.get_cmap('Spectral')
# plot the figure
plt.imshow(ndvi,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
# colour bar
plt.colorbar()

nc.variables.keys()

# better have a look at the individual bands as well
# plot the vis and nir bands
plt.figure(figsize=(8,8))
plt.title('VIS: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(vis,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

plt.figure(figsize=(8,8))
plt.title('NIR: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(nir,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

mask = np.array(nc.variables['Data_Mask']).astype(bool)

# plot it
plt.figure(figsize=(8,8))
plt.title('Data_Mask: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(mask,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

import numpy.ma as ma

vis = ma.array(vis,mask=~mask)
nir = ma.array(nir,mask=~mask)
ndvi = (nir - vis)/(nir + vis)

plt.figure(figsize=(8,8))
plt.title('NDVI: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(ndvi,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

mask1 = vis < 0.
mask2 = nir < 0
print 'number of -ve VIS pixels',np.sum(mask1)
print 'number of -ve NIR pixels',np.sum(mask2)

mask = np.array(nc.variables['Data_Mask']).astype(bool) & (vis > 0) & (nir > 0)

# plot it
plt.figure(figsize=(8,8))
plt.title('Data_Mask: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(mask,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

vis = ma.array(vis,mask=~mask)
nir = ma.array(nir,mask=~mask)
ndvi = (nir - vis)/(nir + vis)

plt.figure(figsize=(8,8))
plt.title('NDVI: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(ndvi,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

# demonstration of multiple subplots

datasets = np.array([['DHR_VIS','DHR_NIR'],\
                     ['DHR_sigmaVIS','DHR_sigmaNIR'],\
                     ['Data_Mask','Weighted_Number_of_Samples']])

# load up all datasets in dict data
data = {}

dlist = datasets.copy().flatten()

for d in dlist:
    data[d] = np.array(nc.variables[d])

mask =  data['Data_Mask'].astype(bool) &( \
        (data['DHR_VIS'] > 0.) | \
        (data['DHR_NIR'] > 0.))
    
s = datasets.shape

# how big for each subplot ?
big = 5

# set the figure size
plt.figure(figsize=(s[1]*big,s[0]*big))

# colorbars for subplots are a bit tricky
# here's one way of sorting this
# using dataset shapes
from matplotlib import gridspec
gs = gridspec.GridSpec(s[0],s[1])

# colour map
cmap = plt.get_cmap('Spectral')


for i,d0 in enumerate(datasets):
    for j,d in enumerate(d0):

        data[d] = ma.array(data[d],mask=~mask)
        
        axes = plt.subplot(gs[i,j])
        axes.set_title(d)
        # no axis ticks
        axes.set_xticks([])
        axes.set_yticks([])
        
        im = axes.imshow(data[d],cmap=cmap,interpolation='none',vmin=0.)
        plt.colorbar(im)

        
ndvi = (data['DHR_NIR'] - data['DHR_VIS'])/(data['DHR_NIR'] + data['DHR_VIS'])

plt.figure(figsize=(8,8))
plt.title('NDVI: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(ndvi,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

# demonstration of multiple subplots

datasets = np.array([['DHR_VIS','DHR_NIR'],\
                     ['DHR_sigmaVIS','DHR_sigmaNIR'],\
                     ['Data_Mask','Weighted_Number_of_Samples']])

# load up all datasets in dict data
data = {}

dlist = datasets.copy().flatten()

for d in dlist:
    data[d] = np.array(nc.variables[d])

mask =  data['Data_Mask'].astype(bool) & \
        (data['Weighted_Number_of_Samples'] > 0.5) & \
        (data['DHR_sigmaVIS'] <= 0.8) & \
        (data['DHR_sigmaNIR'] <= 0.8) & \
        (data['DHR_VIS'] >= 0.) & \
        (data['DHR_NIR'] >= 0.)
    
s = datasets.shape

# how big for each subplot ?
big = 5

# set the figure size
plt.figure(figsize=(s[1]*big,s[0]*big))

# colorbars for subplots are a bit tricky
# here's one way of sorting this
# using dataset shapes
from matplotlib import gridspec
gs = gridspec.GridSpec(s[0],s[1])

# colour map
cmap = plt.get_cmap('Spectral')


for i,d0 in enumerate(datasets):
    for j,d in enumerate(d0):

        data[d] = ma.array(data[d],mask=~mask)
        
        axes = plt.subplot(gs[i,j])
        axes.set_title(d)
        # no axis ticks
        axes.set_xticks([])
        axes.set_yticks([])
        
        im = axes.imshow(data[d],cmap=cmap,interpolation='none',vmin=0.)
        plt.colorbar(im)

ndvi = (data['DHR_NIR'] - data['DHR_VIS'])/(data['DHR_NIR'] + data['DHR_VIS'])

plt.figure(figsize=(13,13))
plt.title('NDVI: Tile %s %d doy %03d'%(tile,year,doy))
cmap = plt.get_cmap('Spectral')
plt.imshow(ndvi,interpolation='none',cmap=cmap,vmin=0.,vmax=1.)
plt.colorbar()

!easy_install --user pysolar

# from https://github.com/pingswept/pysolar/wiki/examples
import Pysolar
from datetime import datetime

# UCL lat/lon
lat = 51.5248
lon = -0.1336

hour = 12
minute = 0
second = 0
month = 10 # ie October
day = 13 
year = 2013

d = datetime(year, month, day, hour, minute, second) 
altitude_deg = Pysolar.GetAltitude(lat, lon, d)
zenith = 90. - altitude_deg
# W m^-2
solar = Pysolar.solar.radiation.GetRadiationDirect(d, altitude_deg)
print zenith,solar

def solar(year, month, day, hour, lat_deg, lon_deg, minute=0, second=0):
    '''Return solar zenith and clear sky radiation
       for given lat, lon and time/date
    '''
    from datetime import datetime
    import Pysolar
    
    d = datetime(year, month, day, hour, minute, second) 
    altitude_deg = Pysolar.GetAltitude(lat_deg, lon_deg, d)
    # W m^-2
    solar_rad = Pysolar.solar.radiation.GetRadiationDirect(d, altitude_deg)
    return 90. - altitude_deg,solar_rad

# or import from local module

import sys,os
# put local directory into the path
sys.path.insert(0,os.path.abspath('files%spython'%os.sep))

from solar import solar

import numpy as np

# UCL lat/lon
lat = 51.5248
lon = -0.1336

second = 0
month = 10 # ie October
day = 13 
year = 2013

radiation_fields = '#hour zenith solar_rad month day lat lon'

radiation = []
for hour in xrange(24):
    for minute in xrange(60):
        thr = hour + minute/60.
        # append data line as tuple
        radiation.append((thr,) + \
                         solar(year, month, day, hour, lat, lon, minute=minute) +\
                         (month, day, lat, lon))
# convert to numpy array
# transpose so access eg zenith as 
# radiation[0]
radiation = np.array(radiation).T

# so we have radiation as 
print radiation.shape
print radiation.ndim
print radiation

import pylab as plt

plt.title('Solar radiation at UCL')
plt.xlabel('hour')
plt.ylabel('solar radiation / Wm^-2')
plt.plot(radiation[0],radiation[2])

import pylab as plt

plt.title('Solar zenith UCL')
plt.xlabel('hour')
plt.ylabel('solar zenith / degrees')
plt.plot(radiation[0],radiation[1])

power_density = radiation[2].sum() * 60
print 'power per unit area = %.3f MJ / m^2'%(power_density/10**6)

import numpy as np

def radiation(year, month, day, lat, lon,minute_step=30):
    rad = []
    for hour in xrange(24):
        for minute in xrange(0,60,minute_step):
            thr = hour + minute/60.
            # append data line as tuple
            rad.append((thr,) + \
                     solar(year, month, day, hour, lat, lon, minute=minute) +\
                     (month, day, lat, lon))
    # convert to numpy array
    # transpose so access eg zenith as 
    # rad[0]
    rad = np.array(rad).T
    return rad

def days_in_month(month,year=2013):
    ''' number of days in month'''
    import calendar
    return calendar.monthrange(year,month)[1]

# UCL lat/lon
lat = 51.5248
lon = -0.1336

year = 2013

minute_step = 30

pd = []
for month in xrange(12):
    ndays = days_in_month(month+1,year=year)
    print month,ndays
    for day in xrange(ndays):
        rad = radiation(year, month+1, day+1, lat, lon,minute_step=minute_step)
        pd.append([month+day/float(ndays),rad[2].sum() * 60 * minute_step])
pd = np.array(pd).T

import pylab as plt

plt.title('Power per unit area, UCL')
plt.xlabel('month')
plt.ylabel('Power per unit area / MJ m^-2')
plt.plot(pd[0],pd[1]/10**6)

# UCL lat/lon
lat = 51.5248
lon = -0.1336

year = 2013

pd = []

for month in xrange(12):
    pd_month = []
    ndays = days_in_month(month+1,year=year)
    print month,ndays
    for day in xrange(ndays):
        rad = radiation(year, month+1, day+1, lat, lon)
        pd_month.append([rad[2].sum() * 60 * 30])
    pd_month = np.array(pd_month).T
    pd.append([month,pd_month.sum()])
pd = np.array(pd).T

import pylab as plt

plt.title('Monthly total Power per unit area, UCL')
plt.xlabel('month')
plt.ylabel('Power per unit area / MJ m^-2')
plt.plot(pd[0],pd[1]/10**6)