|
15 | 15 |
|
16 | 16 | import numpy as np
|
17 | 17 | import pandas as pd
|
| 18 | +import pytz |
18 | 19 |
|
19 | 20 | try:
|
20 | 21 | from .spa_c_files.spa_py import spa_calc
|
|
28 | 29 |
|
29 | 30 |
|
30 | 31 |
|
31 |
| -def get_solarposition(time, location, method='pyephem', pressure=101325, |
| 32 | +def get_solarposition(time, location, method='spa', pressure=101325, |
32 | 33 | temperature=12):
|
33 | 34 | """
|
34 | 35 | A convenience wrapper for the solar position calculators.
|
@@ -199,6 +200,183 @@ def pyephem(time, location, pressure=101325, temperature=12):
|
199 | 200 |
|
200 | 201 |
|
201 | 202 |
|
| 203 | +def ephemeris(time, location, pressure=101325, temperature=12): |
| 204 | + ''' |
| 205 | + Python-native solar position calculator. |
| 206 | + The accuracy of this code is not guaranteed. |
| 207 | + Consider using the built-in spa_c code or the PyEphem library. |
| 208 | +
|
| 209 | + Parameters |
| 210 | + ---------- |
| 211 | + time : pandas.DatetimeIndex |
| 212 | + location : pvlib.Location |
| 213 | +
|
| 214 | + Other Parameters |
| 215 | + ---------------- |
| 216 | +
|
| 217 | + pressure : float or DataFrame |
| 218 | + Ambient pressure (Pascals) |
| 219 | +
|
| 220 | + tempreature: float or DataFrame |
| 221 | + Ambient temperature (C) |
| 222 | + |
| 223 | + Returns |
| 224 | + ------- |
| 225 | + |
| 226 | + DataFrame with the following columns |
| 227 | + |
| 228 | + SunEl : float of DataFrame |
| 229 | + Actual elevation (not accounting for refraction)of the sun |
| 230 | + in decimal degrees, 0 = on horizon. The complement of the True Zenith |
| 231 | + Angle. |
| 232 | + |
| 233 | + SunAz : Azimuth of the sun in decimal degrees from North. 0 = North to 270 = West |
| 234 | + |
| 235 | + SunZen : Solar zenith angle |
| 236 | +
|
| 237 | + ApparentSunEl : float or DataFrame |
| 238 | +
|
| 239 | + Apparent sun elevation accounting for atmospheric |
| 240 | + refraction. This is the complement of the Apparent Zenith Angle. |
| 241 | +
|
| 242 | + SolarTime : fload or DataFrame |
| 243 | + Solar time in decimal hours (solar noon is 12.00). |
| 244 | + |
| 245 | +
|
| 246 | + References |
| 247 | + ----------- |
| 248 | +
|
| 249 | + Grover Hughes' class and related class materials on Engineering |
| 250 | + Astronomy at Sandia National Laboratories, 1985. |
| 251 | +
|
| 252 | + See also |
| 253 | + -------- |
| 254 | + pvl_makelocationstruct |
| 255 | + pvl_alt2pres |
| 256 | + pvl_getaoi |
| 257 | + pvl_spa |
| 258 | +
|
| 259 | + ''' |
| 260 | + |
| 261 | + # Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014 |
| 262 | + # Edited by Will Holmgren (@wholmgren), University of Arizona, 2014 |
| 263 | + |
| 264 | + pvl_logger.warning('Using solarposition.ephemeris is discouraged. ' \ |
| 265 | + 'solarposition.pyephem and solarposition.spa are ' \ |
| 266 | + 'faster and more accurate.') |
| 267 | + |
| 268 | + pvl_logger.debug('location={}, temperature={}, pressure={}'.format(location, temperature, pressure)) |
| 269 | + |
| 270 | + Latitude = location.latitude |
| 271 | + ''' the inversion of longitude is due to the fact that this code was |
| 272 | + originally written for the convention that positive longitude were for |
| 273 | + locations west of the prime meridian. However, the correct convention (as |
| 274 | + of 2009) is to use negative longitudes for locations west of the prime |
| 275 | + meridian. Therefore, the user should input longitude values under the |
| 276 | + correct convention (e.g. Albuquerque is at -106 longitude), but it needs |
| 277 | + to be inverted for use in the code. |
| 278 | + ''' |
| 279 | + Latitude = location.latitude |
| 280 | + Longitude=1 * location.longitude |
| 281 | + Year = time.year |
| 282 | + Month = time.month |
| 283 | + Day = time.day |
| 284 | + Hour = time.hour |
| 285 | + Minute = time.minute |
| 286 | + Second = time.second |
| 287 | + DayOfYear = time.dayofyear |
| 288 | + |
| 289 | + DecHours = Hour + Minute / float(60) + Second / float(3600) |
| 290 | + |
| 291 | + Abber=20 / float(3600) |
| 292 | + LatR=np.radians(Latitude) |
| 293 | + |
| 294 | + # this code is needed to handle the new location.tz format string |
| 295 | + try: |
| 296 | + utc_offset = pytz.timezone(location.tz).utcoffset(time[0]).total_seconds() / 3600. |
| 297 | + except ValueError: |
| 298 | + utc_offset = time.tzinfo.utcoffset(time[0]).total_seconds() / 3600. |
| 299 | + pvl_logger.debug('utc_offset={}'.format(utc_offset)) |
| 300 | + |
| 301 | + UnivDate=DayOfYear |
| 302 | + UnivHr = DecHours + utc_offset - .5 # Will H: surprised to see the 0.5 here, but moving on... |
| 303 | + #+60/float(60)/2 |
| 304 | + |
| 305 | + Yr=Year - 1900 |
| 306 | + |
| 307 | + YrBegin=365 * Yr + np.floor((Yr - 1) / float(4)) - 0.5 |
| 308 | + |
| 309 | + Ezero=YrBegin + UnivDate |
| 310 | + T=Ezero / float(36525) |
| 311 | + GMST0=6 / float(24) + 38 / float(1440) + (45.836 + 8640184.542 * T + 0.0929 * T ** 2) / float(86400) |
| 312 | + GMST0=360 * (GMST0 - np.floor(GMST0)) |
| 313 | + GMSTi=np.mod(GMST0 + 360 * (1.0027379093 * UnivHr / float(24)),360) |
| 314 | + |
| 315 | + LocAST=np.mod((360 + GMSTi - Longitude),360) |
| 316 | + EpochDate=Ezero + UnivHr / float(24) |
| 317 | + T1=EpochDate / float(36525) |
| 318 | + ObliquityR=np.radians(23.452294 - 0.0130125 * T1 - 1.64e-06 * T1 ** 2 + 5.03e-07 * T1 ** 3) |
| 319 | + MlPerigee=281.22083 + 4.70684e-05 * EpochDate + 0.000453 * T1 ** 2 + 3e-06 * T1 ** 3 |
| 320 | + MeanAnom=np.mod((358.47583 + 0.985600267 * EpochDate - 0.00015 * T1 ** 2 - 3e-06 * T1 ** 3),360) |
| 321 | + Eccen=0.01675104 - 4.18e-05 * T1 - 1.26e-07 * T1 ** 2 |
| 322 | + EccenAnom=MeanAnom |
| 323 | + E=0 |
| 324 | + |
| 325 | + while np.max(abs(EccenAnom - E)) > 0.0001: |
| 326 | + E=EccenAnom |
| 327 | + EccenAnom=MeanAnom + np.degrees(Eccen)*(np.sin(np.radians(E))) |
| 328 | + |
| 329 | + #pdb.set_trace() |
| 330 | + TrueAnom=2 * np.mod(np.degrees(np.arctan2(((1 + Eccen) / (1 - Eccen)) ** 0.5*np.tan(np.radians(EccenAnom) / float(2)),1)),360) |
| 331 | + EcLon=np.mod(MlPerigee + TrueAnom,360) - Abber |
| 332 | + EcLonR=np.radians(EcLon) |
| 333 | + DecR=np.arcsin(np.sin(ObliquityR)*(np.sin(EcLonR))) |
| 334 | + Dec=np.degrees(DecR) |
| 335 | + #pdb.set_trace() |
| 336 | + |
| 337 | + RtAscen=np.degrees(np.arctan2(np.cos(ObliquityR)*((np.sin(EcLonR))),np.cos(EcLonR))) |
| 338 | + |
| 339 | + HrAngle=LocAST - RtAscen |
| 340 | + HrAngleR=np.radians(HrAngle) |
| 341 | + |
| 342 | + HrAngle=HrAngle - (360*((abs(HrAngle) > 180))) |
| 343 | + SunAz=np.degrees(np.arctan2(- 1 * np.sin(HrAngleR),np.cos(LatR)*(np.tan(DecR)) - np.sin(LatR)*(np.cos(HrAngleR)))) |
| 344 | + SunAz=SunAz + (SunAz < 0) * 360 |
| 345 | + SunEl=np.degrees(np.arcsin((np.cos(LatR)*(np.cos(DecR))*(np.cos(HrAngleR)) + np.sin(LatR)*(np.sin(DecR))))) #potential error |
| 346 | + SolarTime=(180 + HrAngle) / float(15) |
| 347 | + |
| 348 | + Refract=[] |
| 349 | + |
| 350 | + for Elevation in SunEl: |
| 351 | + TanEl=np.tan(np.radians(Elevation)) |
| 352 | + if Elevation>5 and Elevation<=85: |
| 353 | + Refract.append((58.1 / float(TanEl) - 0.07 / float(TanEl ** 3) + 8.6e-05 / float(TanEl ** 5))) |
| 354 | + elif Elevation > -0.575 and Elevation <=5: |
| 355 | + Refract.append((Elevation*((- 518.2 + Elevation*((103.4 + Elevation*((- 12.79 + Elevation*(0.711))))))) + 1735)) |
| 356 | + elif Elevation> -1 and Elevation<= -0.575: |
| 357 | + Refract.append(- 20.774 / float(TanEl)) |
| 358 | + else: |
| 359 | + Refract.append(0) |
| 360 | + |
| 361 | + |
| 362 | + Refract=np.array(Refract)*((283 / float(273 + temperature)))*(pressure) / float(101325) / float(3600) |
| 363 | + |
| 364 | + SunZen = 90 - SunEl |
| 365 | + #SunZen[SunZen >= 90] = 90 |
| 366 | + |
| 367 | + ApparentSunEl = SunEl + Refract |
| 368 | + |
| 369 | + DFOut = pd.DataFrame({'elevation':SunEl}, index=time) |
| 370 | + DFOut['azimuth'] = SunAz-180 #Changed RA Feb 18,2014 to match Duffe |
| 371 | + DFOut['zenith'] = SunZen |
| 372 | + DFOut['apparent_elevation'] = ApparentSunEl |
| 373 | + DFOut['apparent_zenith'] = 90 - ApparentSunEl |
| 374 | + DFOut['solar_time'] = SolarTime |
| 375 | + |
| 376 | + return DFOut |
| 377 | + |
| 378 | + |
| 379 | + |
202 | 380 | def _localize_to_utc(time, location):
|
203 | 381 | """
|
204 | 382 | Converts or localizes a time series to UTC.
|
|
0 commit comments