/
julian.py
325 lines (259 loc) · 8.15 KB
/
julian.py
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
"""
Created on Tue Apr 02 16:50:34 2013
@author: tm
computes julian date, given month (1..12). day(1..31) and year,
its inverse (calendar date from julian), and the day
of the year (doy), assuming it is a leap year.
julday and caldat are adapted from "Numerical Recipes in C', 2nd edition,
pp. 11
restrictions
- no error handling implemented
- works only for years past 1582
- time not yet supported
"""
import numpy as np
import pandas as pd
import datetime as dt
import pytz
def julday(month, day, year, hour=0, minute=0, second=0):
"""
Julian date from month, day, and year (can be scalars or arrays)
Parameters
----------
month : numpy.ndarray or int32
Month.
day : numpy.ndarray or int32
Day.
year : numpy.ndarray or int32
Year.
hour : numpy.ndarray or int32, optional
Hour.
minute : numpy.ndarray or int32, optional
Minute.
second : numpy.ndarray or int32, optional
Second.
Returns
-------
jul : numpy.ndarray or double
Julian day.
"""
month = np.array(month)
day = np.array(day)
inJanFeb = month <= 2
jy = year - inJanFeb
jm = month + 1 + inJanFeb * 12
jul = np.int32(np.floor(365.25 * jy) +
np.floor(30.6001 * jm) + (day + 1720995.0))
ja = np.int32(0.01 * jy)
jul += 2 - ja + np.int32(0.25 * ja)
jul = jul + hour / 24.0 - 0.5 + minute / 1440.0 + second / 86400.0
return jul
def caldat(julian):
"""
Calendar date (month, day, year) from julian date, inverse of 'julday()'
Return value: month, day, and year in the Gregorian
Works only for years past 1582!
Parameters
----------
julian : numpy.ndarray or double
Julian day.
Returns
-------
month : numpy.ndarray or int32
Month.
day : numpy.ndarray or int32
Day.
year : numpy.ndarray or int32
Year.
"""
jn = np.int32(((np.array(julian) + 0.000001).round()))
jalpha = np.int32(((jn - 1867216) - 0.25) / 36524.25)
ja = jn + 1 + jalpha - (np.int32(0.25 * jalpha))
jb = ja + 1524
jc = np.int32(6680.0 + ((jb - 2439870.0) - 122.1) / 365.25)
jd = np.int32(365.0 * jc + (0.25 * jc))
je = np.int32((jb - jd) / 30.6001)
day = jb - jd - np.int32(30.6001 * je)
month = je - 1
month = (month - 1) % 12 + 1
year = jc - 4715
year = year - (month > 2)
return month, day, year
def julian2date(julian):
"""
Calendar date from julian date.
Works only for years past 1582!
Parameters
----------
julian : numpy.ndarray or double
Julian day.
Returns
-------
year : numpy.ndarray or int32
Year.
month : numpy.ndarray or int32
Month.
day : numpy.ndarray or int32
Day.
hour : numpy.ndarray or int32
Hour.
minute : numpy.ndarray or int32
Minute.
second : numpy.ndarray or int32
Second.
"""
min_julian = 2299160
max_julian = 1827933925
julian = np.atleast_1d(np.array(julian, dtype=float))
if np.min(julian) < min_julian or np.max(julian) > max_julian:
raise ValueError("Value of Julian date is out of allowed range.")
jn = (np.round(julian + 0.0000001)).astype(np.int32)
jalpha = (((jn - 1867216) - 0.25) / 36524.25).astype(np.int32)
ja = jn + 1 + jalpha - (np.int32(0.25 * jalpha))
jb = ja + 1524
jc = (6680.0 + ((jb - 2439870.0) - 122.1) / 365.25).astype(np.int32)
jd = (365.0 * jc + (0.25 * jc)).astype(np.int32)
je = ((jb - jd) / 30.6001).astype(np.int32)
day = jb - jd - np.int64(30.6001 * je)
month = je - 1
month = (month - 1) % 12 + 1
year = jc - 4715
year = year - (month > 2)
fraction = (julian + 0.5 - jn).astype(np.float64)
eps = (np.float64(1e-12) * np.abs(jn)).astype(np.float64)
eps.clip(min=np.float64(1e-12), max=None)
hour = (fraction * 24. + eps).astype(np.int64)
hour.clip(min=0, max=23)
fraction -= hour / 24.
minute = (fraction * 1440. + eps).astype(np.int64)
minute.clip(min=0, max=59)
second = (fraction - minute / 1440.) * 86400.
second.clip(min=0, max=None)
microsecond = ((second - np.int32(second)) * 1e6).astype(np.int32)
second = second.astype(np.int32)
return year, month, day, hour, minute, second, microsecond
def julian2datetime(julian, tz=None):
"""
converts julian date to python datetime
default is not time zone aware
Parameters
----------
julian : float
julian date
"""
year, month, day, hour, minute, second, microsecond = julian2date(julian)
if type(julian) == np.array or type(julian) == np.memmap or \
type(julian) == np.ndarray or type(julian) == np.flatiter:
return np.array([dt.datetime(y, m, d, h, mi, s, ms, tz)
for y, m, d, h, mi, s, ms in
zip(np.atleast_1d(year),
np.atleast_1d(month),
np.atleast_1d(day),
np.atleast_1d(hour),
np.atleast_1d(minute),
np.atleast_1d(second),
np.atleast_1d(microsecond))])
return dt.datetime(year, month, day, hour, minute, second, microsecond, tz)
def julian2doy(j, consider_nonleap_years=True):
"""
Calendar date from julian date.
Works only for years past 1582!
Parameters
----------
j : numpy.ndarray or double
Julian days.
consider_nonleap_years : boolean, optional
Flag if all dates are interpreted as leap years (False) or not (True).
Returns
-------
doy : numpy.ndarray or int32
Day of year.
"""
year, month, day = julian2date(j)[0:3]
if consider_nonleap_years:
return doy(month, day, year)
else:
return doy(month, day)
def julian2datetimeindex(j, tz=pytz.UTC):
"""
Converting Julian days to datetimeindex.
Parameters
----------
j : numpy.ndarray or int32
Julian days.
tz : instance of pytz, optional
Time zone. Default: UTC
Returns
-------
datetime : pandas.DatetimeIndex
Datetime index.
"""
year, month, day, hour, minute, second, microsecond = julian2date(j)
return pd.DatetimeIndex([dt.datetime(y, m, d, h, mi, s, ms, tz)
for y, m, d, h, mi, s, ms in
zip(year, month, day, hour, minute,
second, microsecond)])
def julian2num(j):
"""
Convert a matplotlib date to a Julian days.
Parameters
----------
j : numpy.ndarray : int32
Julian days.
Returns
-------
num : numpy.ndarray : int32
Number of days since 0001-01-01 00:00:00 UTC *plus* *one*.
"""
return j - 1721424.5
def num2julian(n):
"""
Convert a Julian days to a matplotlib date.
Parameters
----------
n : numpy.ndarray : int32
Number of days since 0001-01-01 00:00:00 UTC *plus* *one*.
Returns
-------
j : numpy.ndarray : int32
Julian days.
"""
return n + 1721424.5
def doy(month, day, year=None):
"""
Calculation of day of year. If year is provided it will be tested for
leap years.
Parameters
----------
month : numpy.ndarray or int32
Month.
day : numpy.ndarray or int32
Day.
year : numpy.ndarray or int32, optional
Year.
Retruns
-------
doy : numpy.ndarray or int32
Day of year.
"""
daysPast = np.array([0, 31, 60, 91, 121, 152, 182, 213,
244, 274, 305, 335, 366])
day_of_year = daysPast[month - 1] + day
if year is not None:
nonleap_years = np.invert(is_leap_year(year))
day_of_year = day_of_year - nonleap_years + \
np.logical_and(day_of_year < 60, nonleap_years)
return day_of_year
def is_leap_year(year):
"""
Check if year is a leap year.
Parameters
----------
year : numpy.ndarray or int32
Returns
-------
leap_year : numpy.ndarray or boolean
True if year is a leap year.
"""
return np.logical_or(np.logical_and(year % 4 == 0, year % 100 != 0),
year % 400 == 0)