weixin_39717865
weixin_39717865
2020-12-27 05:17

Preserve seconds during back and forth conversion of Julian dates

Dear Jeff,

I am using the routines date2num and num2date and realised that they do not give the same results when I convert back and forth: In [1]: import netCDF4 as nt In [2]: jd0 = 1721424. In [3]: jd = nt.date2num(nt.datetime(1810, 4, 24, 16, 15, 10), 'days since 0001-01-01 12:00:00', calendar='gregorian')+jd0 In [4]: nt.num2date(jd-jd0, 'days since 0001-01-01 12:00:00', calendar='gregorian') Out[4]: datetime.datetime(1810, 4, 24, 16, 15, 9, 999999)

jd0 is the Julian Day of 0001-01-01 12:00:00. As you can see, I get only 9 seconds back instead of second 10.

In my Fortran code, I normally add a little epsilon so that the back conversion works. When I do this: In [1]: import netCDF4 as nt In [2]: jd0 = 1721424. In [3]: jd = nt.date2num(nt.datetime(1810, 4, 24, 16, 15, 10, 5), 'days since 0001-01-01 12:00:00', calendar='gregorian')+jd0 In [4]: nt.num2date(jd-jd0, 'days since 0001-01-01 12:00:00:00', calendar='gregorian') Out[4]: datetime.datetime(1810, 4, 24, 16, 15, 10)

As you can see, it works well for everything but the microseconds. This is because 1 microsecond is In [5]: 1./1e6/86400. Out[5]: 1.1574074074074074e-11 whereas the correction is ca. In [7]: 2382262.1771990741*2.22E-016 Out[7]: 5.288622033381944e-10 that means around 45 microseconds.

Until now, I do 1. Add the epsilon (above) in JulianDayFromDate. 2. Do not add 0.00005 to Z in DateFromJulianDay als set eps = 0. 3. Substract the epsilon only from the microseconds not on anything before.

With this I am always less then 30 microseconds high. With your formulation, we are also less than 30 microseconds off but with a tendency for lower numbers. This is undesireable because 0 microseconds is the default (and what most people would want) and your code gives often the second just below the real second.

My code is always a bit high so that the seconds flips erronously at about 999980 microseconds.

If one wants more precision, one has to go to 128-bit. One could borrow from the astropy library, which does exactly this and gets nanosecond resolution.

Please find my changes. Until now there is a global switch domc = True that switches between old (False) and new (True) formulations.

Kind regards Matthias

该提问来源于开源项目:Unidata/netcdf4-python

  • 点赞
  • 写回答
  • 关注问题
  • 收藏
  • 复制链接分享
  • 邀请回答

6条回答

  • weixin_39658900 weixin_39658900 4月前

    one of the tests in tst_netcdftime.py is failing....

    点赞 评论 复制链接分享
  • weixin_39717865 weixin_39717865 4月前

    Dear Jeff,

    this is because t = 733498.999999 is 8640 microseconds off the real 733499. You applied the numerical fix in the back conversion to real dates while I apply it during the conversion to the Julian date. So this test will not work anymore. However, adding 4 more 9s to the end fixes it: t = 733498.9999999999 Now the test succeeds because then it is less then 1 microsecond off the real 733499.

    But I have to admit that now another test fails: # issue 415 t = datetimex(2001, 12, 1, 2, 3, 4) self.assertEqual(t, copy.deepcopy(t)) I would claim that I have nothing to do with it ;-) Must have been failing before.

    Regards Matthias

    On 26 Jun 2015, at 15:06, Jeff Whitaker notifications.com wrote:

    one of the tests in tst_netcdftime.py is failing....

    — Reply to this email directly or view it on GitHub.

    点赞 评论 复制链接分享
  • weixin_39658900 weixin_39658900 4月前

    Looks like all the tests are passing now, thanks.

    点赞 评论 复制链接分享
  • weixin_39658900 weixin_39658900 4月前

    I think this can be merged after the change suggested above (using np.finfo instead of hard-coded eps value). I also suggest removing the if domc blocks - this makes the code harder to read and it's easy enough to revert the change later if we need to.

    点赞 评论 复制链接分享
  • weixin_39658900 weixin_39658900 4月前

    Also, please add a Changelog entry.

    点赞 评论 复制链接分享
  • weixin_39658900 weixin_39658900 4月前

    Thanks - merging now...

    点赞 评论 复制链接分享

相关推荐