-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcode
123 lines (112 loc) · 4.59 KB
/
code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
# -*- coding: utf-8 -*-
# Solar zenith angle calculation model
# author:shaodonghang
# date:2017/12/4
# Reference:王炳忠.太阳辐射计算讲座:第一讲 太阳能中天文参数的计算[J].太阳能,1999(2):8-10.
import time
import numpy as np
import datetime as datetime
import matplotlib.pylab as plt
import netCDF4 as nc
import sys
import math
import xlrd
# import gdal
from xlwt import *
from xlutils.copy import copy
# import gdal
import ogr
from gdalconst import *
from osr import SpatialReference
from osgeo import gdal
from osgeo.gdalconst import *
# import libtiff as tif
# from libtiff import TIFF #read tiff data
from datetime import timedelta
from scipy.interpolate import griddata
from netCDF4 import Dataset
#########################################
data = xlrd.open_workbook('C:/Users/sdh/Desktop/SolarAngle.xlsx')
# read time and lon/lat data
table = data.sheets()[0] # 打开第一张表
nrows = table.nrows # 获取表的行数
ncols = table.ncols
station = table.col_values(0)[1:]
year = table.col_values(1)[1:]
month = table.col_values(2)[1:]
day = table.col_values(3)[1:]
hour = table.col_values(4)[1:]
minute = table.col_values(5)[1:]
sec = table.col_values(6)[1:]
lon = table.col_values(7)[1:]
lat = table.col_values(8)[1:]
TimeZone = table.col_values(9)[1:]
#################################
solarAngle=[]
wb = copy(data)
ws = wb.get_sheet(0)
ws.write(0, 11, 'Day of Year')
ws.write(0, 12, 'Local Time') #实际是gtdt
ws.write(0, 13, 'Sun Angle')#(sitar)
ws.write(0, 14, 'Declination Angle')
ws.write(0, 15, 'Equation of Time')
style = easyxf('pattern: pattern solid, fore_color yellow;')
ws.write(0, 16, 'ZenithAngle(deg)',style)
ws.write(0, 17, 'HeightAngle(deg)',style)
ws.write(0, 18, 'AzimuthAngle(deg)',style)
for n in range(1,nrows):
m=n-1
#年积日的计算
#儒略日 Julian day(由通用时转换到儒略日)
JD0 = int(365.25*(year[m]-1))+int(30.6001*(1+13))+1+hour[m]/24+1720981.5
if month[m]<=2:
JD2 = int(365.25*(year[m]-1))+int(30.6001*(month[m]+13))+day[m]+hour[m]/24+1720981.5
else:
JD2 = int(365.25*year[m])+int(30.6001*(month[m]+1))+day[m]+hour[m]/24+1720981.5
#年积日 Day of year
DOY = JD2-JD0+1
#N0 sitar=θ
N0 = 79.6764 + 0.2422*(year[m]-1985) - int((year[m]-1985)/4.0)
sitar = 2*math.pi*(DOY-N0)/365.2422
ED1 = 0.3723 + 23.2567*math.sin(sitar) + 0.1149*math.sin(2*sitar) - 0.1712*math.sin(3*sitar)- 0.758*math.cos(sitar) + 0.3656*math.cos(2*sitar) + 0.0201*math.cos(3*sitar)
ED = ED1*math.pi/180 #ED本身有符号
if lon[m] >= 0:
if TimeZone == -13:
dLon = lon[m] - (math.floor((lon[m]*10-75)/150)+1)*15.0
else:
dLon = lon[m] - TimeZone[m]*15.0 #地球上某一点与其所在时区中心的经度差
else:
if TimeZone[m] == -13:
dLon = (math.floor((lon[m]*10-75)/150)+1)*15.0- lon[m]
else:
dLon = TimeZone[m]*15.0- lon[m]
#时差
Et = 0.0028 - 1.9857*math.sin(sitar) + 9.9059*math.sin(2*sitar) - 7.0924*math.cos(sitar)- 0.6882*math.cos(2*sitar)
gtdt1 = hour[m] + minute[m]/60.0 + sec[m]/3600.0 + dLon/15 #地方时
gtdt = gtdt1 + Et/60.0
dTimeAngle1 = 15.0*(gtdt-12)
dTimeAngle = dTimeAngle1*math.pi/180
latitudeArc = lat[m]*math.pi/180
# 高度角计算公式
HeightAngleArc = math.asin(math.sin(latitudeArc)*math.sin(ED)+math.cos(latitudeArc)*math.cos(ED)*math.cos(dTimeAngle))
# 方位角计算公式
CosAzimuthAngle = (math.sin(HeightAngleArc)*math.sin(latitudeArc)-math.sin(ED))/math.cos(HeightAngleArc)/math.cos(latitudeArc)
AzimuthAngleArc = math.acos(CosAzimuthAngle)
HeightAngle = HeightAngleArc*180/math.pi
ZenithAngle = 90-HeightAngle
AzimuthAngle1 = AzimuthAngleArc *180/math.pi
# print 'ZenithAngle=',ZenithAngle
# print 'HeightAngle=',HeightAngle
if dTimeAngle < 0:
AzimuthAngle = 180 - AzimuthAngle1
else:
AzimuthAngle = 180 + AzimuthAngle1
# print 'AzimuthAngle=',AzimuthAngle
print('station:'+station[m]+' '+'ZenithAngle(deg):%f HeightAngle(deg):%f AzimuthAngle(deg):%f' % (ZenithAngle,HeightAngle,AzimuthAngle))
# np.savetxt(outpath +'CL_UHB' +str(i)+'.txt',inew_CL_nc )
solarAngle.append(ZenithAngle)
# solarAngle=np.array(solarAngle)
np.savetxt('C:/Users/sdh/Desktop/result/SolarAngle.txt',solarAngle)
# print 'ZenithAngle=',ZenithAngle
# print 'HeightAngle=',HeightAngle
# print 'AzimuthAngle=',AzimuthAngle