def absorbed_power_density(year,minute_step=60): '''Function to calculate the monthly total absorbed solar radiation power density in MJ / m^2 for a given year given an input dataset albedo and associated lat and lon information (in degrees). The shape of the albedo dataset is: (12,nlat,nlon) Aguments: year : integer of the year Options: minute_step : integer: resolution of steps in minutes. Must be a divisor of 60 (e.g. 10, 15, 30, 60) ''' import sys sys.path.insert(0,'files/python') # modification of masked.py from masked2 import masked import numpy as np year = 2010 data = masked(dataset=['lat','lon','BHR_SW'],year=year) # take the first lat and lon only as they are all the same lat = data['lat'][0] lon = data['lon'][0] albedo = data['BHR_SW'] s = albedo.shape print s # pad out the lat and lon arrays # this is a bit tricky so test this in parts # we dont need this repeated for each month though lat2 = np.array([lat] * s[2]).T lon2 = np.array([lon] * s[1]) print lat2.shape print lon2.shape import Pysolar from datetime import datetime 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 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 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 # functions so far import sys sys.path.insert(0,'files/python') # modification of masked.py from masked2 import masked import numpy as np import Pysolar from datetime import datetime 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 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] ''' Processing -- test for several lat/long ''' import numpy.ma as ma year = 2010 minute_step = 60 data = masked(dataset=['lat','lon','BHR_SW'],year=year) lat = data['lat'][0] lon = data['lon'][0] albedo = data['BHR_SW'] s = albedo.shape lat = ma.array( [lat] *s[2]).T lon = ma.array([lon] * s[1]) # array for output absorbed = ma.array(np.zeros_like(albedo),mask=albedo.mask) # #for i in xrange(lat.shape[0]): # for j in xrange(lat.shape[1]): # try it out for 6 samples in longitude, including first and last lons = range(0,lon.shape[1],lon.shape[1]/5) + [lon.shape[1]-1] # an array to hold the sample pd for the sample # longitudes all_pd = [] i = lon.shape[0]/2 for j in lons: pd = [] # loop over month for month in xrange(12): pd_month = [] ndays = days_in_month(month+1,year=year) print i,j,month,ndays # loop over days for day in xrange(ndays): # solar radiation for that day rad = radiation(year, month+1, day+1, \ lat[i,j], lon[i,j],minute_step=minute_step) pd_month.append([rad[2].sum() * 60 * minute_step]) pd_month = np.array(pd_month).T pd.append([pd_month.sum()]) # pd is the power density for each month all_pd.append(pd) # MJ per m2 all_pd = np.array(all_pd).squeeze() import pylab as plt plt.plot(lon[0,lons],all_pd) plt.xlim(-180.,180.) plt.xlabel('latitude') plt.ylabel('power density total per month') def absorbed_power_density(year,minute_step=60): '''Function to calculate the monthly total absorbed solar radiation power density in MJ / m^2 for a given year given an input dataset albedo and associated lat and lon information (in degrees). The shape of the albedo dataset is: (12,nlat,nlon) Aguments: year : integer of the year Options: minute_step : integer: resolution of steps in minutes. Must be a divisor of 60 (e.g. 10, 15, 30, 60) ''' data = masked(dataset=['lat','lon','BHR_SW'],year=year) lat = data['lat'][0] lon = data['lon'][0] albedo = data['BHR_SW'] s = albedo.shape lat = ma.array( [lat] *s[2]).T lon = ma.array([lon] * s[1]) # array for output absorbed = ma.array(np.zeros_like(albedo),mask=albedo.mask) power_density = ma.array(np.zeros_like(albedo),mask=albedo.mask) # for i in xrange(0,lat.shape[0]): print 'lat',lat[i,0] # single longitude for j in [lat.shape[1]/2]: pd = [] # loop over month for month in xrange(12): pd_month = [] ndays = days_in_month(month+1,year=year) # loop over days for day in xrange(ndays): # solar radiation for that day rad = radiation(year, month+1, day+1, \ lat[i,j], lon[i,j],minute_step=minute_step) # get rid of -ves rad[rad<0] = 0 pd_month.append([rad[2].sum() * 60 * minute_step]) pd_month = np.array(pd_month).T pd.append([pd_month.sum()]) pd = np.array(pd).squeeze().T # pd is the power density for each month # load this into all longitudes for k in xrange(power_density.shape[0]): power_density[k,i,:] = pd[k] power_density /= 10.**6 # MJ per m2 absorbed = power_density * (1 - albedo) return absorbed,power_density # this will still take some minutes to calculate absorbed,power_density = absorbed_power_density(year,minute_step=60) # this is how you can save a numpy array np.savez('files/data/absorbed.npz',absorbed=np.array(absorbed),\ power_density=np.array(power_density),mask=absorbed.mask) # now next time we should be able to load it try: f = np.load('files/data/absorbed.npz') mask = f['mask'] absorbed,power_density = ma.array(f['absorbed'],mask=mask),\ ma.array(f['power_density'],mask=mask) except: absorbed,power_density = absorbed_power_density(year,minute_step=60) import pylab as plt cmap = plt.get_cmap('Spectral') vmax = absorbed.max()/2 vmin = 0. for m in xrange(12): plt.figure(figsize=(7,3)) plt.title('Absorbed radiation Month %02d %d MJ / m^2'%(m,year)) plt.imshow(absorbed[m],cmap=cmap,interpolation='none',vmin=vmin,vmax=vmax) plt.colorbar() plt.savefig('files/data/absorbed%02d.jpg'%m) import pylab as plt cmap = plt.get_cmap('Spectral') vmax = absorbed.max() * 2./3. vmin = 0. for m in xrange(12): plt.figure(figsize=(7,3)) plt.title('Solar radiation power density Month %02d %d MJ / m^2'%(m,year)) plt.imshow(power_density[m],cmap=cmap,interpolation='none',vmin=vmin,vmax=vmax) plt.colorbar() plt.savefig('files/data/power%02d.jpg'%m) # make some movies!??? import os for t in ['power','absorbed']: for m in xrange(12): cmd = 'convert files/data/%s%02d.jpg files/data/%s%02d.gif'%(t,m,t,m) print cmd os.system(cmd) cmd = "convert -delay 100 -loop 0 \ files/data/%s??.gif files/data/%s_movie.gif"%(t,t) os.system(cmd)