ls -l data/*.nc

import glob

files = glob.glob('data/*.nc')
print files

filename = 'data/data.dat'

# open in write mode
fp = open(filename,'w')
fp.writelines(files)
fp.close()

cat data/data.dat

files = glob.glob('data/*.nc')

for i,file in enumerate(files):
    files[i] = file + '\n'
print files

files = glob.glob('data/*.nc')

# or:
files = [file + '\n' for file in files]

print files

# or all at once if you like:

files = [file + '\n' for file in glob.glob('data/*.nc')]

print files

import glob

files = [file + '\n' for file in glob.glob('data/*.nc')]

filename = 'data/data.dat'

# open in write mode
fp = open(filename,'w')
fp.writelines(files)
fp.close()

cat data/data.dat

import gdal

# form a generic string of the form
# NETCDF:"data/GlobAlbedo.200901.mosaic.5.nc":DHR_VIS

file_template = 'NETCDF:"%s":%s'

# now make a list of the datset names we want
# so we can loop over this 

selected_layers = ['DHR_VIS','DHR_NIR','DHR_SW']

# ----------------------------------

# try it out:


root = 'data/'

# example filename : use formatting string:
# %d%02d
year = 2009
month = 1
filename = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)
print filename


# set up an empty dictionary to load the data into
data = {}

# use enumerate here to loop over
# the list selected_layers and also have
# access to an index i

for i, layer in enumerate ( selected_layers ):
    this_file = file_template % ( filename, layer )
    print "Opening Layer %d: %s" % (i, this_file )
    g = gdal.Open ( this_file )
    
    # test that the opening worked
    # and raise an error otherwise
    
    if g is None:
        raise IOError
    data[layer] = g.ReadAsArray() 
    print "\t>>> Read %s!" % layer


# post-hoc load into array called albedo
albedo = [data[f] for f in selected_layers]


import gdal

# NETCDF:"data/GlobAlbedo.200901.mosaic.5.nc":DHR_VIS
file_template = 'NETCDF:"%s":%s'

selected_layers = ['DHR_VIS','DHR_NIR','DHR_SW']

root = 'data/'

# example filename : use formatting string:
# %d%02d
year = 2009
month = 1
filename = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)


# set up an empty list to load the data into
data = []

# loop over
# the list selected_layers 

for layer in ( selected_layers ):
    this_file = file_template % ( filename, layer )
    g = gdal.Open ( this_file )
    
    # test that the opening worked
    # and raise an error otherwise
    
    if g is None:
        raise IOError
        
    # it is a list, so we will just append the entry here
    data.append(g.ReadAsArray()) 
    
# check it workwed
print 'previously, we had:',albedo
print '...... we now have:',data

import gdal

# NETCDF:"data/GlobAlbedo.200901.mosaic.5.nc":DHR_VIS
file_template = 'NETCDF:"%s":%s'

selected_layers = ['DHR_VIS','DHR_NIR','DHR_SW']

root = 'data/'

# example filename : use formatting string:
# %d%02d
year = 2009
month = 1
filename = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)


# set up an empty list to load the data into
data = []

# loop over
# the list selected_layers 

for layer in ( selected_layers ):
    this_file = file_template % ( filename, layer )
    g = gdal.Open ( this_file )
    
    # test that the opening worked
    # and raise an error otherwise
    
    if g is None:
        raise IOError
        
    # it is a list, so we will just append the entry here
    data.append(g.ReadAsArray()) 

import gdal

# NETCDF:"data/GlobAlbedo.200901.mosaic.5.nc":DHR_VIS
file_template = 'NETCDF:"%s":%s'

selected_layers = ['DHR_VIS','DHR_NIR','DHR_SW']

root = 'data/'

# example filename : use formatting string:
# %d%02d
year = 2009

# important to put the alldata empty list
# setup outside of the loop, ie before
# we start looping over month

alldata = []

for month in range(1,13):
  print "I'm reading month",month,"of year",year
  filename = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)
   
  # just the data for this month in the list data  
  data = []

  for layer in ( selected_layers ):
    this_file = file_template % ( filename, layer )
    g = gdal.Open ( this_file )
    
    # test that the opening worked
    # and raise an error otherwise
    
    if g is None:
        raise IOError
        
    # it is a list, so we will just append the entry here
    data.append(g.ReadAsArray()) 
  # now (note indentation!!)
  # append to the alldata list
  alldata.append(data)
    
# now check how big ...

print 'dataset dimensions',len(alldata),'x',len(alldata[0]),'x',\
                len(alldata[0][0]),'x',len(alldata[0][0][0])

!convert data/albedo.jpg data/albedo.gif

import os
cmd = 'convert data/albedo.jpg data/albedo.gif'
os.system(cmd)

import gdal
import pylab as plt
import os
import calendar

layer = 'BHR_VIS'

# form a generic string of the form
# NETCDF:"data/GlobAlbedo.200901.mosaic.5.nc":DHR_VIS

file_template = 'NETCDF:"%s":%s'

root = 'data/'

# example filename : use formatting string:
# %d%02d
year = 2009

# set the month (1-based, i.e. 1 == January)
month = 1

filename = 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()


''' Plot the data and save as picture jpeg format '''
# make a string with the output file name
out_file = root + 'GlobAlbedo.%d%02d.jpg'%(year,month)
# plot
plt.figure(figsize=(10,4))
plt.clf()
# %9s forces the string to be 8 characters long
plt.title('VIS BHR albedo for %8s %d'%(calendar.month_name[month],year))
# use nearest neighbour interpolation
# load the array data 
plt.imshow(data,interpolation='nearest',cmap=plt.get_cmap('Spectral'),vmin=0.0,vmax=1.0)
# show a colour bar 
plt.colorbar()
plt.savefig(out_file)

# convert to gif
# set up the unix command which is of the form
# convert input output
# Here input will be out_file
# and output we can get with out_file.replace('.jpg','.gif')
# i.e. replacing where it says .jpg with .gif
cmd = 'convert %s %s'%(out_file,out_file.replace('.jpg','.gif'))
os.system(cmd)


ls -l data/GlobAlbedo.??????.gif

# lets try something out to see how we can loop easily in this case
# we know that the mionth names are contained in calendar.month_name
# so we might try to just loop over that

import calendar

for month_name in calendar.month_name:
  print month_name

# but we might also want access to the month index, so we might use
# enumerate

import calendar

for month,month_name in enumerate(calendar.month_name):
  print month,month_name

# but thats not quite right ... we don't want
# the blank zero entry so we only
# loop over the slice [1:] in month_name

import calendar

for month,month_name in enumerate(calendar.month_name[1:]):
  print month,month_name

# but now we have a zero index for January
# which isnt what we want
# we could just add one to this when we use it
# but that is a bit ugly ...
# so instead, use a feature of enumerate()

import calendar

for month,month_name in enumerate(calendar.month_name[1:],start=1):
  print month,month_name

import gdal
import pylab as plt
import os
import calendar

# set the layer we want
layer = 'BHR_SW'

# form a generic string of the form
# NETCDF:"data/GlobAlbedo.200901.mosaic.5.nc":BHR_SW

file_template = 'NETCDF:"%s":%s'

root = 'data/'

# example filename : use formatting string:
# %d%02d
year = 2009

# set the month (1-based, i.e. 1 == January)
# its a good idea to use enumerate here 
# but we want the numbers to start at 1
# so we set start=1 in the enumerate call
for month,month_name in enumerate(calendar.month_name[1:],start=1):
  '''
  Loop over each month
  
  setting month = 1,2,..12
          month_name= 'January', ... 'December'
          (unless we change the locale!)
  '''
    
  filename = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)
  g = gdal.Open (  file_template % ( filename, layer ) )

  if g is None:
   raise IOError
   
  # note that we are re-using the varioable data
  # each time in the loop as we haver no need to store it as yet
        
  data = g.ReadAsArray()

  ''' Plot the data and save as picture jpeg format '''
  # make a string with the output file name
  out_file = root + 'GlobAlbedo.%d%02d.jpg'%(year,month)
  # plot  
  plt.figure(figsize=(10,4))
  plt.clf()
  # %9s forces the string to be 8 characters long
  plt.title('SW BHR albedo for %8s %d'%(month_name,year))
  # use nearest neighbour interpolation
  # load the array data 
  plt.imshow(data,interpolation='nearest',cmap=plt.get_cmap('Spectral'),vmin=0.0,vmax=1.0)
  # show a colour bar 
  plt.colorbar()
  plt.savefig(out_file)

  # convert to gif
  # set up the unix command which is of the form
  # convert input output
  # Here input will be out_file
  # and output we can get with out_file.replace('.jpg','.gif')
  # i.e. replacing where it says .jpg with .gif
  cmd = 'convert %s %s'%(out_file,out_file.replace('.jpg','.gif')) 
  os.system(cmd)


ls -l data/GlobAlbedo.2009.SW.gif

import os

# this is quite a neat way of generating the string for the input files
out_file = 'data/GlobAlbedo.%d.SW.gif'%year
in_files = out_file.replace('.SW.gif','??.gif')

cmd = 'convert -delay 100 -loop 0 %s %s'%(in_files,out_file)
# check the cmd is ok
print cmd

# good ... so now run it
os.system(cmd)

from netCDF4 import Dataset
import numpy as np

root = 'files/data/'
year = 2009

# which months?
months = xrange(1,13)

# empty list
data = []

# loop over month
# use enumerate so we have an index counter
for i,month in enumerate(months):
    # this then is the file we want
    local_file = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)
    
    # load the netCDF data from the file local_file
    nc = Dataset(local_file,'r')
    # append what we read to the list called data
    data.append(np.array(nc.variables['DHR_SW']))
    
# convert data to a numpy array (its a list of arrays at the moment)
data = np.array(data)

import numpy.ma as ma

band = np.array(nc.variables['DHR_SW'])

masked_band = ma.array(band,mask=np.isnan(band))
print masked_band.mask

from netCDF4 import Dataset
import numpy as np
import numpy.ma as ma

root = 'files/data/'
year = 2009

# which months?
months = xrange(1,13)

# empty list
data = []

# loop over month
# use enumerate so we have an index counter
for i,month in enumerate(months):
    # this then is the file we want
    local_file = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)
    
    # load the netCDF data from the file local_file
    nc = Dataset(local_file,'r')
    # load into the variable 'band'
    band = np.array(nc.variables['DHR_SW'])
    # convert to a masked array
    masked_band = ma.array(band,mask=np.isnan(band))
    # append what we read to the list called data
    data.append(masked_band)
    
# convert data to a numpy array (its a list of arrays at the moment)
data = np.array(data)

print type(data)
print data.shape
print data.ndim

from netCDF4 import Dataset
import numpy as np
import numpy.ma as ma

root = 'files/data/'
year = 2009

# which months?
months = xrange(1,13)

# empty list
data = []

# loop over month
# use enumerate so we have an index counter
for i,month in enumerate(months):
    # this then is the file we want
    local_file = root + 'GlobAlbedo.%d%02d.mosaic.5.nc'%(year,month)
    
    # load the netCDF data from the file local_file
    nc = Dataset(local_file,'r')
    # load into the variable 'band'
    band = np.array(nc.variables['DHR_SW'])
    # convert to a masked array
    masked_band = ma.array(band,mask=np.isnan(band))
    # append what we read to the list called data
    data.append(masked_band)
    
# convert data to a numpy array (its a list of arrays at the moment)
data = ma.array(data)

print type(data)
print data.shape
print data.ndim

print data.mask.shape
print data.mask.ndim

sum = (data.mask).sum(axis=0)
print sum.shape

print sum

np.unique(sum)