准备数据
# 经度、维度、时区
longitude = 115.81567872975157
dimensionality = 28.900022800650184
timezone = 8
# 以下是js用于计算日出日落时间
mjd_def = """
function mjd(day, month, year, hour) {
var a, b;
if (month <= 2) {
month = month + 12;
year = year - 1;
}
a = 10000.0 * year + 100.0 * month + day;
if (a <= 15821004.1) {
b = -2 * Math.floor((year + 4716)/4) - 1179;
}
else {
b = Math.floor(year/400) - Math.floor(year/100) + Math.floor(year/4);
}
a = 365.0 * year - 679004.0;
return (a + b + Math.floor(30.6001 * (month + 1)) + day + hour/24.0);
}
"""
cal_def = """
function hrsmin(hours) {
var hrs, h, m, dum;
hrs = Math.floor(hours * 60 + 0.5)/ 60.0;
h = Math.floor(hrs);
m = Math.floor(60 * (hrs - h) + 0.5);
if(h<10)h="0"+h;
if(m<10)m="0"+m;
dum = h +":"+ m;
if (dum < 1000) dum = "0" + dum;
if (dum <100) dum = "0" + dum;
if (dum < 10) dum = "0" + dum;
return dum;
}
function ipart(x) {
var a;
if (x> 0) {
a = Math.floor(x);
}
else {
a = Math.ceil(x);
}
return a;
}
function frac(x) {
var a;
a = x - Math.floor(x);
if (a < 0) a += 1;
return a;
}
function round(num, dp) {
return Math.round (num * Math.pow(10, dp)) / Math.pow(10, dp);
}
function range(x) {
var a, b;
b = x / 360;
a = 360 * (b - ipart(b));
if (a < 0 ) {
a = a + 360
}
return a
}
function mjd(day, month, year, hour) {
var a, b;
if (month <= 2) {
month = month + 12;
year = year - 1;
}
a = 10000.0 * year + 100.0 * month + day;
if (a <= 15821004.1) {
b = -2 * Math.floor((year + 4716)/4) - 1179;
}
else {
b = Math.floor(year/400) - Math.floor(year/100) + Math.floor(year/4);
}
a = 365.0 * year - 679004.0;
return (a + b + Math.floor(30.6001 * (month + 1)) + day + hour/24.0);
}
function quad(ym, yz, yp) {
var nz, a, b, c, dis, dx, xe, ye, z1, z2, nz;
var quadout = new Array;
nz = 0;
a = 0.5 * (ym + yp) - yz;
b = 0.5 * (yp - ym);
c = yz;
xe = -b / (2 * a);
ye = (a * xe + b) * xe + c;
dis = b * b - 4.0 * a * c;
if (dis > 0) {
dx = 0.5 * Math.sqrt(dis) / Math.abs(a);
z1 = xe - dx;
z2 = xe + dx;
if (Math.abs(z1) <= 1.0) nz += 1;
if (Math.abs(z2) <= 1.0) nz += 1;
if (z1 < -1.0) z1 = z2;
}
quadout[0] = nz;
quadout[1] = z1;
quadout[2] = z2;
quadout[3] = xe;
quadout[4] = ye;
return quadout;
}
function lmst(mjday, glong) {
var lst, t, d;
d = mjday - 51544.5
t = d / 36525.0;
lst = range(280.46061837 + 360.98564736629 * d + 0.000387933 *t*t - t*t*t / 38710000);
return (lst/15.0 + glong/15);
}
function minisun(t) {
var p2 = 6.283185307, coseps = 0.91748, sineps = 0.39778;
var L, M, DL, SL, X, Y, Z, RHO, ra, dec;
var suneq = new Array;
M = p2 * frac(0.993133 + 99.997361 * t);
DL = 6893.0 * Math.sin(M) + 72.0 * Math.sin(2 * M);
L = p2 * frac(0.7859453 + M / p2 + (6191.2 * t + DL)/1296000);
SL = Math.sin(L);
X = Math.cos(L);
Y = coseps * SL;
Z = sineps * SL;
RHO = Math.sqrt(1 - Z * Z);
dec = (360.0 / p2) * Math.atan(Z / RHO);
ra = (48.0 / p2) * Math.atan(Y / (X + RHO));
if (ra <0 ) ra += 24;
suneq[1] = dec;
suneq[2] = ra;
return suneq;
}
function sin_alt(iobj, mjd0, hour, glong, cglat, sglat) {
var mjday, t, ra, dec, tau, salt, rads = 0.0174532925;
var objpos = new Array;
mjday = mjd0 + hour/24.0;
t = (mjday - 51544.5) / 36525.0;
if (iobj == 1) {
objpos = minimoon(t);
}
else {
objpos = minisun(t);
}
ra = objpos[2];
dec = objpos[1];
tau = 15.0 * (lmst(mjday, glong) - ra);
salt = sglat * Math.sin(rads*dec) + cglat * Math.cos(rads*dec) * Math.cos(rads*tau);
return salt;
}
function getzttime(mjday, tz, glong) {
var sglong, sglat, date, ym, yz, utrise, utset, j;
var yp, nz, hour, z1, z2, iobj, rads = 0.0174532925;
var quadout = new Array;
sinho = Math.sin(rads * -0.833);
date = mjday - tz/24;
hour = 1.0;
ym = sin_alt(2, date, hour - 1.0, glong, 1, 0) - sinho;
while(hour < 25) {
yz = sin_alt(2, date, hour, glong, 1, 0) - sinho;
yp = sin_alt(2, date, hour + 1.0, glong, 1, 0) - sinho;
quadout = quad(ym, yz, yp);
nz = quadout[0];
z1 = quadout[1];
z2 = quadout[2];
xe = quadout[3];
ye = quadout[4];
if (nz == 1) {
if (ym < 0.0)
utrise = hour + z1;
else
utset = hour + z1;
}
if (nz == 2) {
if (ye < 0.0) {
utrise = hour + z2;
utset = hour + z1;
}
else {
utrise = hour + z1;
utset = hour + z2;
}
}
ym = yp;
hour += 2.0;
}
var zt=(utrise+utset)/2;
if(zt<utrise)
zt=(zt+12)%24;
return zt;
}
function Cal(mjday, tz, glong, glat) {
var sglong, sglat, date, ym, yz, above, utrise, utset, j;
var yp, nz, rise, sett, hour, z1, z2, iobj, rads = 0.0174532925;
var quadout = new Array;
var always_up = "日不落";
var always_down = "日不出";
var outstring = "";
var resobj = {};
sinho = Math.sin(rads * -0.833);
sglat = Math.sin(rads * glat);
cglat = Math.cos(rads * glat);
date = mjday - tz/24;
rise = false;
sett = false;
above = false;
hour = 1.0;
ym = sin_alt(2, date, hour - 1.0, glong, cglat, sglat) - sinho;
if (ym > 0.0) above = true;
while(hour < 25 && (sett == false || rise == false)) {
yz = sin_alt(2, date, hour, glong, cglat, sglat) - sinho;
yp = sin_alt(2, date, hour + 1.0, glong, cglat, sglat) - sinho;
quadout = quad(ym, yz, yp);
nz = quadout[0];
z1 = quadout[1];
z2 = quadout[2];
xe = quadout[3];
ye = quadout[4];
if (nz == 1) {
if (ym < 0.0) {
utrise = hour + z1;
rise = true;
}
else {
utset = hour + z1;
sett = true;
}
}
if (nz == 2) {
if (ye < 0.0) {
utrise = hour + z2;
utset = hour + z1;
}
else {
utrise = hour + z1;
utset = hour + z2;
}
}
ym = yp;
hour += 2.0;
}
if (rise == true || sett == true ) {
if (rise == true) {
resobj["rise"] = hrsmin(utrise);
outstring += " 日出时刻:" + hrsmin(utrise)+"<BR>";
}
else{
outstring += " 日不出或日不落"+"<BR>";
resobj["pole"] = "日不出或日不落";
}
var zt=getzttime(mjday, tz, glong);
resobj["center"] = hrsmin(zt);
outstring+= " 日中时刻:" + hrsmin(zt)+"<BR>";
if (sett == true){
outstring += " 日没时刻:" + hrsmin(utset);
resobj["set"] = hrsmin(utset);
}else {
outstring += " 日不出或日不落";
resobj["pole"] = "日不出或日不落";
}
}
else {
if (above == true){
outstring += always_up;
var zt=getzttime(mjday, tz, glong);
resobj["center"] = hrsmin(zt);
resobj["pole"] = "极昼";
outstring+="<BR>"+" 日中时刻:"+hrsmin(zt);
} else {
outstring += always_down;
resobj["pole"] = "极夜";
}
}
return resobj;
}
"""
调用js计算
import execjs
def calculate_rise_set_center(datatime,timezone,longitude_latitude):
"""
parameter:
----------
datatime: tuple(day, month, year, hour)
timezone: int 8
longitude_latitude: tuple(longitude, latitude)
====================
return:
dict(rise, center, set)
example:
calculate_rise_set_center((12, 10, 2018, 0.0),8,(116.3833, 39.9))
return: {'rise': '06:21', 'center': '12:01', 'set': '17:41'}
"""
# 计算mjd
ctx = execjs.compile(mjd_def)
mjd = ctx.call("mjd", datatime[0], datatime[1], datatime[2], datatime[3]) # 传入：天、月、年、小时
# 计算日出日落
ctx = execjs.compile(cal_def)
result = ctx.call("Cal", mjd, timezone, longitude_latitude[0], longitude_latitude[1]) # 传入：时区、经度、维度
return result
import datetime
day = datetime.datetime.now().day
month = datetime.datetime.now().month
year = datetime.datetime.now().year
print(calculate_rise_set_center((day, month, year, 0.0),timezone,(longitude, dimensionality)))
{'rise': '06:18', 'center': '12:02', 'set': '17:46'}
计算结果
print(day_list)
print(rise_list)
print(set_list)
print(center_list)
[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]
['06:09', '06:10', '06:10', '06:11', '06:12', '06:12', '06:13', '06:13', '06:14', '06:14', '06:15', '06:15', '06:16', '06:17', '06:17', '06:18', '06:18', '06:19', '06:20', '06:20', '06:21', '06:22', '06:22', '06:23', '06:24', '06:24', '06:25', '06:26', '06:26', '06:27', '06:28']
['18:03', '18:02', '18:01', '18:00', '17:58', '17:57', '17:56', '17:55', '17:54', '17:53', '17:52', '17:51', '17:50', '17:49', '17:48', '17:46', '17:45', '17:44', '17:43', '17:42', '17:41', '17:41', '17:40', '17:39', '17:38', '17:37', '17:36', '17:35', '17:34', '17:33', '17:33']
['12:06', '12:06', '12:06', '12:05', '12:05', '12:05', '12:05', '12:04', '12:04', '12:04', '12:04', '12:03', '12:03', '12:03', '12:03', '12:02', '12:02', '12:02', '12:02', '12:02', '12:01', '12:01', '12:01', '12:01', '12:01', '12:01', '12:01', '12:01', '12:00', '12:00', '12:00']
绘图
from matplotlib import pyplot
# 显示中文,ubuntu mono字体
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
font_20 = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=20)
font_50 = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=50)
# 设置图片大小
pyplot.figure(figsize=(20, 8), dpi=100)
# 画图，day_list做X轴，rise_list做Y轴, 这条线表示日出, 用#FF3E96表示, 实线，实心转折点，4倍的线宽， 50%的透明度
pyplot.plot(day_list, rise_list, label='日出', color="#FF3E96", linestyle="solid", marker='o', linewidth=4, alpha=0.5)
# 画第二条线，这条线表示日落
pyplot.plot(day_list, set_list, label='日落', marker='x')
#标题
pyplot.title("2018年10月日出时间", fontproperties=font_20)
# X轴表示天数，Y轴表示时间
pyplot.xlabel("日期", fontproperties=font)
pyplot.ylabel("时间", fontproperties=font)
# X轴刻度按day_list显示，但每3个显示一个，可读性高，Y轴也一样
pyplot.xticks(day_list[::2])
# 格式化Y轴
# rise_show = ["{}点{}分".format(rise_time.split(":")[0], rise_time.split(":")[1]) for rise_time in rise_list]
# pyplot.yticks(rise_list[::4], rise_show[::4], fontproperties=font, rotation=10) # 添加字体, 旋转
# 合并两个Y轴
rise_set_list = rise_list + set_list
rise_set_show = ["{}点{}分".format(rise_set_time.split(":")[0], rise_set_time.split(":")[1]) for rise_set_time in rise_set_list]
pyplot.yticks(rise_set_list[::4], rise_set_show[::4], fontproperties=font)
# 显示网格,默认每个刻度,alpha透明度,网格线
pyplot.grid(alpha=0.5, linestyle=':')
# 显示图例,设置字体,位置
pyplot.legend(prop=font, loc="lower right")
# 标识最大值
# 求最大索引
max_index = rise_list.index(max(rise_list))
# 在图上画点
pyplot.plot(max_index+1,rise_list[max_index],'ks')
# 标记的文本
show_max='['+str(max_index+1)+' '+str(rise_list[max_index])+']'
# 把文本标记上去， xytext标记在哪个坐标， xy表示标记哪个点
pyplot.annotate(show_max,xytext=(max_index-1,rise_list[max_index]),xy=(max_index,rise_list[max_index]))
# 加水印
pyplot.text(1, "17:37", r'Walker出品',alpha=0.5,fontproperties=font_50, color='red')
# 保存图片
pyplot.savefig('./lesson_1_1.png')
walker日出日落图
