@@ -23,6 +23,10 @@ from PIL import Image from math import * import numpy as np +import solar +import datetime as dt +import simulate as sim +import radiation as rad def squareImage(im): (width, height) = im.size @@ -40,13 +44,28 @@ def despherifyImage(im): full_circle = 1000.0 * 2 * pi +# for r in range(half_width): +# for theta in range(int(full_circle)): +# (inx, iny) = (round(r * cos(theta/1000.0)) + half_width, round(r * sin(theta/1000.0)) + half_width) +# (outx, outy) = (width - width * (theta/full_circle) - 1, r) +# outpix[outx, outy] = inpix[inx, iny] + + theta_range = range(int(full_circle)) + t_1000_range = [t / 1000.0 for t in theta_range] + thetas = zip([cos(t) for t in t_1000_range], + [sin(t) for t in t_1000_range], + [t / full_circle for t in theta_range]) + for r in range(half_width): - for theta in range(int(full_circle)): - (inx, iny) = (round(r * cos(theta/1000.0)) + half_width, round(r * sin(theta/1000.0)) + half_width) - (outx, outy) = (width - width * (theta/full_circle) - 1, r) - outpix[outx, outy] = inpix[inx, iny] + for t_cos, t_sin, t_full_circle in thetas: + (inx, iny) = (round(r * t_cos) + half_width, + round(r * t_sin) + half_width) + outx = width - width * (t_full_circle) - 1 + outpix[outx, r] = inpix[inx, iny] return out + + def differentiateImageColumns(im): (width, height) = im.size pix = im.load() @@ -62,21 +81,69 @@ def redlineImage(im): pix = im.load() threshold = 300 + horizon = [] for x in range(width): for y in range(height - 1): (R, G, B) = pix[x, y] if R + G + B > threshold: pix[x, y] = (255, 0, 0) + horizon.append(y) break + return (im, horizon) + +def addSunPaths(im, latitude_deg, longitude_deg, horizon, d): + pix = im.load() + + alt_zero = getAltitudeZero() + az_zero = getAzimuthZero() + + for m in range(24 * 60): + alt = solar.GetAltitude(latitude_deg, longitude_deg, d + dt.timedelta(minutes = m)) + az = solar.GetAzimuth(latitude_deg, longitude_deg, d + dt.timedelta(minutes = m)) + + x = az_zero + int(az * 1944.0/360.0) + y = alt_zero - int(alt_zero * sin(radians(alt))) + if y < horizon[x]: + pix[x % 1944, y] = (255, 193, 37) + return im +def getAzimuthZero(): + return 1100 + +def getAltitudeZero(): + return 380 + +horizon = [] + im = Image.open('spherical.jpg').convert("L") im = squareImage(im) +print 'Starting despherification . . .' lin = despherifyImage(im) + +print 'Despherification complete. Calculating horizon . . .' d = differentiateImageColumns(lin).convert("RGB") -r = redlineImage(d) +r, horizon = redlineImage(d) +print 'Horizon calculated.' + +(latitude_deg, longitude_deg) = (42.206, -71.382) +summer = dt.datetime(2009, 6, 21, 5, 0, 0, 0) +fall = dt.datetime(2009, 9, 21, 5, 0, 0, 0) +winter = dt.datetime(2009, 12, 21, 5, 0, 0, 0) +step_minutes = 5 + +power_densities = [radiation for (time, alt, az, radiation, shade) in sim.SimulateSpan(latitude_deg, longitude_deg, horizon, summer, winter, step_minutes)] +print power_densities + +energy = sum(power_densities) * step_minutes * 60 +print str(energy/1000000) + ' MJ per m^2 per year' + +sp = addSunPaths(r, latitude_deg, longitude_deg, horizon, summer) +sp2 = addSunPaths(r, latitude_deg, longitude_deg, horizon, fall) +sp3 = addSunPaths(r, latitude_deg, longitude_deg, horizon, winter) -r.show() +sp3.show() +#sp3.save('sun_path.jpg') horizon.py @@ -22,6 +22,7 @@ import datetime import radiation import solar +from math import * def BuildTimeList(start_utc_datetime, end_utc_datetime, step_minutes): '''Create a list of sample points evenly spaced apart by step_minutes.''' @@ -29,10 +30,18 @@ def BuildTimeList(start_utc_datetime, end_utc_datetime, step_minutes): time_list = [] span = end_utc_datetime - start_utc_datetime dt = datetime.timedelta(seconds = step) - print span return map(lambda n: start_utc_datetime + dt * n, range((span.days * 86400 + span.seconds) / step)) -def SimulateSpan(latitude_deg, longitude_deg, start_utc_datetime, end_utc_datetime, step_minutes, elevation = 0, temperature_celsius = 25, pressure_millibars = 1013.25): +def CheckAgainstHorizon(power): + (time, alt, az, radiation, shade) = power + alt_zero = 380 + + if shade < alt_zero - int(alt_zero * sin(radians(alt))): + radiation = 0 + + return (time, alt, az, radiation, shade) + +def SimulateSpan(latitude_deg, longitude_deg, horizon, start_utc_datetime, end_utc_datetime, step_minutes, elevation = 0, temperature_celsius = 25, pressure_millibars = 1013.25): '''Simulate the motion of the sun over a time span and location of your choosing. The start and end points are set by datetime objects, which can be created with @@ -47,8 +56,8 @@ def SimulateSpan(latitude_deg, longitude_deg, start_utc_datetime, end_utc_dateti solar.GetAltitude(latitude_deg, longitude_deg, time, elevation, temperature_celsius, pressure_millibars), solar.GetAzimuth(latitude_deg, longitude_deg, time, elevation) ) for time in time_list] - power_list = [(time, alt, az, radiation.GetRadiationDirect(time, alt)) for (time, alt, az) in angles_list] - print power_list + power_list = [(time, alt, az, radiation.GetRadiationDirect(time, alt), horizon[int(az)]) for (time, alt, az) in angles_list] + return filter(CheckAgainstHorizon, power_list) # xs = shade.GetXShade(width, 120, azimuth_deg) # ys = shade.GetYShade(height, 120, altitude_deg) simulate.py 6b8bc63da99523ae2dd892d98354e19d6f42a6cb Brandon Stafford brandon@pingswept.org http://github.com/pingswept/pysolar/commit/2d5b33c7c4561950e9895d072f40543034389772 2d5b33c7c4561950e9895d072f40543034389772 2009-04-10T20:44:58-07:00 2009-04-10T20:44:58-07:00 Added code for overlaying sun paths. 7c9a18ea1a3f2e24395b4d2d2d5ad67984ddc1e9 Brandon Stafford brandon@pingswept.org