1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2""" 3The astropy.utils.iers package provides access to the tables provided by 4the International Earth Rotation and Reference Systems Service, in 5particular allowing interpolation of published UT1-UTC values for given 6times. These are used in `astropy.time` to provide UT1 values. The polar 7motions are also used for determining earth orientation for 8celestial-to-terrestrial coordinate transformations 9(in `astropy.coordinates`). 10""" 11 12import re 13from datetime import datetime 14from warnings import warn 15from urllib.parse import urlparse 16 17import numpy as np 18import erfa 19 20from astropy.time import Time, TimeDelta 21from astropy import config as _config 22from astropy import units as u 23from astropy.table import QTable, MaskedColumn 24from astropy.utils.data import (get_pkg_data_filename, clear_download_cache, 25 is_url_in_cache, get_readable_fileobj) 26from astropy.utils.state import ScienceState 27from astropy import utils 28from astropy.utils.exceptions import AstropyWarning 29 30__all__ = ['Conf', 'conf', 'earth_orientation_table', 31 'IERS', 'IERS_B', 'IERS_A', 'IERS_Auto', 32 'FROM_IERS_B', 'FROM_IERS_A', 'FROM_IERS_A_PREDICTION', 33 'TIME_BEFORE_IERS_RANGE', 'TIME_BEYOND_IERS_RANGE', 34 'IERS_A_FILE', 'IERS_A_URL', 'IERS_A_URL_MIRROR', 'IERS_A_README', 35 'IERS_B_FILE', 'IERS_B_URL', 'IERS_B_README', 36 'IERSRangeError', 'IERSStaleWarning', 37 'LeapSeconds', 'IERS_LEAP_SECOND_FILE', 'IERS_LEAP_SECOND_URL', 38 'IETF_LEAP_SECOND_URL'] 39 40# IERS-A default file name, URL, and ReadMe with content description 41IERS_A_FILE = 'finals2000A.all' 42IERS_A_URL = 'ftp://anonymous:mail%40astropy.org@gdc.cddis.eosdis.nasa.gov/pub/products/iers/finals2000A.all' # noqa: E501 43IERS_A_URL_MIRROR = 'https://datacenter.iers.org/data/9/finals2000A.all' 44IERS_A_README = get_pkg_data_filename('data/ReadMe.finals2000A') 45 46# IERS-B default file name, URL, and ReadMe with content description 47IERS_B_FILE = get_pkg_data_filename('data/eopc04_IAU2000.62-now') 48IERS_B_URL = 'http://hpiers.obspm.fr/iers/eop/eopc04/eopc04_IAU2000.62-now' 49IERS_B_README = get_pkg_data_filename('data/ReadMe.eopc04_IAU2000') 50 51# LEAP SECONDS default file name, URL, and alternative format/URL 52IERS_LEAP_SECOND_FILE = get_pkg_data_filename('data/Leap_Second.dat') 53IERS_LEAP_SECOND_URL = 'https://hpiers.obspm.fr/iers/bul/bulc/Leap_Second.dat' 54IETF_LEAP_SECOND_URL = 'https://www.ietf.org/timezones/data/leap-seconds.list' 55 56# Status/source values returned by IERS.ut1_utc 57FROM_IERS_B = 0 58FROM_IERS_A = 1 59FROM_IERS_A_PREDICTION = 2 60TIME_BEFORE_IERS_RANGE = -1 61TIME_BEYOND_IERS_RANGE = -2 62 63MJD_ZERO = 2400000.5 64 65INTERPOLATE_ERROR = """\ 66interpolating from IERS_Auto using predictive values that are more 67than {0} days old. 68 69Normally you should not see this error because this class 70automatically downloads the latest IERS-A table. Perhaps you are 71offline? If you understand what you are doing then this error can be 72suppressed by setting the auto_max_age configuration variable to 73``None``: 74 75 from astropy.utils.iers import conf 76 conf.auto_max_age = None 77""" 78 79MONTH_ABBR = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 80 'Sep', 'Oct', 'Nov', 'Dec'] 81 82 83def download_file(*args, **kwargs): 84 """ 85 Overload astropy.utils.data.download_file within iers module to use a 86 custom (longer) wait time. This just passes through ``*args`` and 87 ``**kwargs`` after temporarily setting the download_file remote timeout to 88 the local ``iers.conf.remote_timeout`` value. 89 """ 90 kwargs.setdefault('http_headers', {'User-Agent': 'astropy/iers', 91 'Accept': '*/*'}) 92 93 with utils.data.conf.set_temp('remote_timeout', conf.remote_timeout): 94 return utils.data.download_file(*args, **kwargs) 95 96 97def _none_to_float(value): 98 """ 99 Convert None to a valid floating point value. Especially 100 for auto_max_age = None. 101 """ 102 return (value if value is not None else np.finfo(float).max) 103 104 105class IERSStaleWarning(AstropyWarning): 106 pass 107 108 109class Conf(_config.ConfigNamespace): 110 """ 111 Configuration parameters for `astropy.utils.iers`. 112 """ 113 auto_download = _config.ConfigItem( 114 True, 115 'Enable auto-downloading of the latest IERS data. If set to False ' 116 'then the local IERS-B file will be used by default (even if the ' 117 'full IERS file with predictions was already downloaded and cached). ' 118 'This parameter also controls whether internet resources will be ' 119 'queried to update the leap second table if the installed version is ' 120 'out of date. Default is True.') 121 auto_max_age = _config.ConfigItem( 122 30.0, 123 'Maximum age (days) of predictive data before auto-downloading. ' 124 'See "Auto refresh behavior" in astropy.utils.iers documentation for details. ' 125 'Default is 30.') 126 iers_auto_url = _config.ConfigItem( 127 IERS_A_URL, 128 'URL for auto-downloading IERS file data.') 129 iers_auto_url_mirror = _config.ConfigItem( 130 IERS_A_URL_MIRROR, 131 'Mirror URL for auto-downloading IERS file data.') 132 remote_timeout = _config.ConfigItem( 133 10.0, 134 'Remote timeout downloading IERS file data (seconds).') 135 system_leap_second_file = _config.ConfigItem( 136 '', 137 'System file with leap seconds.') 138 iers_leap_second_auto_url = _config.ConfigItem( 139 IERS_LEAP_SECOND_URL, 140 'URL for auto-downloading leap seconds.') 141 ietf_leap_second_auto_url = _config.ConfigItem( 142 IETF_LEAP_SECOND_URL, 143 'Alternate URL for auto-downloading leap seconds.') 144 145 146conf = Conf() 147 148 149class IERSRangeError(IndexError): 150 """ 151 Any error for when dates are outside of the valid range for IERS 152 """ 153 154 155class IERS(QTable): 156 """Generic IERS table class, defining interpolation functions. 157 158 Sub-classed from `astropy.table.QTable`. The table should hold columns 159 'MJD', 'UT1_UTC', 'dX_2000A'/'dY_2000A', and 'PM_x'/'PM_y'. 160 """ 161 162 iers_table = None 163 """Cached table, returned if ``open`` is called without arguments.""" 164 165 @classmethod 166 def open(cls, file=None, cache=False, **kwargs): 167 """Open an IERS table, reading it from a file if not loaded before. 168 169 Parameters 170 ---------- 171 file : str or None 172 full local or network path to the ascii file holding IERS data, 173 for passing on to the ``read`` class methods (further optional 174 arguments that are available for some IERS subclasses can be added). 175 If None, use the default location from the ``read`` class method. 176 cache : bool 177 Whether to use cache. Defaults to False, since IERS files 178 are regularly updated. 179 180 Returns 181 ------- 182 IERS 183 An IERS table class instance 184 185 Notes 186 ----- 187 On the first call in a session, the table will be memoized (in the 188 ``iers_table`` class attribute), and further calls to ``open`` will 189 return this stored table if ``file=None`` (the default). 190 191 If a table needs to be re-read from disk, pass on an explicit file 192 location or use the (sub-class) close method and re-open. 193 194 If the location is a network location it is first downloaded via 195 download_file. 196 197 For the IERS class itself, an IERS_B sub-class instance is opened. 198 199 """ 200 if file is not None or cls.iers_table is None: 201 if file is not None: 202 if urlparse(file).netloc: 203 kwargs.update(file=download_file(file, cache=cache)) 204 else: 205 kwargs.update(file=file) 206 207 # TODO: the below is really ugly and probably a bad idea. Instead, 208 # there should probably be an IERSBase class, which provides 209 # useful methods but cannot really be used on its own, and then 210 # *perhaps* an IERS class which provides best defaults. But for 211 # backwards compatibility, we use the IERS_B reader for IERS here. 212 if cls is IERS: 213 cls.iers_table = IERS_B.read(**kwargs) 214 else: 215 cls.iers_table = cls.read(**kwargs) 216 return cls.iers_table 217 218 @classmethod 219 def close(cls): 220 """Remove the IERS table from the class. 221 222 This allows the table to be re-read from disk during one's session 223 (e.g., if one finds it is out of date and has updated the file). 224 """ 225 cls.iers_table = None 226 227 def mjd_utc(self, jd1, jd2=0.): 228 """Turn a time to MJD, returning integer and fractional parts. 229 230 Parameters 231 ---------- 232 jd1 : float, array, or `~astropy.time.Time` 233 first part of two-part JD, or Time object 234 jd2 : float or array, optional 235 second part of two-part JD. 236 Default is 0., ignored if jd1 is `~astropy.time.Time`. 237 238 Returns 239 ------- 240 mjd : float or array 241 integer part of MJD 242 utc : float or array 243 fractional part of MJD 244 """ 245 try: # see if this is a Time object 246 jd1, jd2 = jd1.utc.jd1, jd1.utc.jd2 247 except Exception: 248 pass 249 250 mjd = np.floor(jd1 - MJD_ZERO + jd2) 251 utc = jd1 - (MJD_ZERO+mjd) + jd2 252 return mjd, utc 253 254 def ut1_utc(self, jd1, jd2=0., return_status=False): 255 """Interpolate UT1-UTC corrections in IERS Table for given dates. 256 257 Parameters 258 ---------- 259 jd1 : float, array of float, or `~astropy.time.Time` object 260 first part of two-part JD, or Time object 261 jd2 : float or float array, optional 262 second part of two-part JD. 263 Default is 0., ignored if jd1 is `~astropy.time.Time`. 264 return_status : bool 265 Whether to return status values. If False (default), 266 raise ``IERSRangeError`` if any time is out of the range covered 267 by the IERS table. 268 269 Returns 270 ------- 271 ut1_utc : float or float array 272 UT1-UTC, interpolated in IERS Table 273 status : int or int array 274 Status values (if ``return_status``=``True``):: 275 ``iers.FROM_IERS_B`` 276 ``iers.FROM_IERS_A`` 277 ``iers.FROM_IERS_A_PREDICTION`` 278 ``iers.TIME_BEFORE_IERS_RANGE`` 279 ``iers.TIME_BEYOND_IERS_RANGE`` 280 """ 281 return self._interpolate(jd1, jd2, ['UT1_UTC'], 282 self.ut1_utc_source if return_status else None) 283 284 def dcip_xy(self, jd1, jd2=0., return_status=False): 285 """Interpolate CIP corrections in IERS Table for given dates. 286 287 Parameters 288 ---------- 289 jd1 : float, array of float, or `~astropy.time.Time` object 290 first part of two-part JD, or Time object 291 jd2 : float or float array, optional 292 second part of two-part JD (default 0., ignored if jd1 is Time) 293 return_status : bool 294 Whether to return status values. If False (default), 295 raise ``IERSRangeError`` if any time is out of the range covered 296 by the IERS table. 297 298 Returns 299 ------- 300 D_x : `~astropy.units.Quantity` ['angle'] 301 x component of CIP correction for the requested times. 302 D_y : `~astropy.units.Quantity` ['angle'] 303 y component of CIP correction for the requested times 304 status : int or int array 305 Status values (if ``return_status``=``True``):: 306 ``iers.FROM_IERS_B`` 307 ``iers.FROM_IERS_A`` 308 ``iers.FROM_IERS_A_PREDICTION`` 309 ``iers.TIME_BEFORE_IERS_RANGE`` 310 ``iers.TIME_BEYOND_IERS_RANGE`` 311 """ 312 return self._interpolate(jd1, jd2, ['dX_2000A', 'dY_2000A'], 313 self.dcip_source if return_status else None) 314 315 def pm_xy(self, jd1, jd2=0., return_status=False): 316 """Interpolate polar motions from IERS Table for given dates. 317 318 Parameters 319 ---------- 320 jd1 : float, array of float, or `~astropy.time.Time` object 321 first part of two-part JD, or Time object 322 jd2 : float or float array, optional 323 second part of two-part JD. 324 Default is 0., ignored if jd1 is `~astropy.time.Time`. 325 return_status : bool 326 Whether to return status values. If False (default), 327 raise ``IERSRangeError`` if any time is out of the range covered 328 by the IERS table. 329 330 Returns 331 ------- 332 PM_x : `~astropy.units.Quantity` ['angle'] 333 x component of polar motion for the requested times. 334 PM_y : `~astropy.units.Quantity` ['angle'] 335 y component of polar motion for the requested times. 336 status : int or int array 337 Status values (if ``return_status``=``True``):: 338 ``iers.FROM_IERS_B`` 339 ``iers.FROM_IERS_A`` 340 ``iers.FROM_IERS_A_PREDICTION`` 341 ``iers.TIME_BEFORE_IERS_RANGE`` 342 ``iers.TIME_BEYOND_IERS_RANGE`` 343 """ 344 return self._interpolate(jd1, jd2, ['PM_x', 'PM_y'], 345 self.pm_source if return_status else None) 346 347 def _check_interpolate_indices(self, indices_orig, indices_clipped, max_input_mjd): 348 """ 349 Check that the indices from interpolation match those after clipping 350 to the valid table range. This method gets overridden in the IERS_Auto 351 class because it has different requirements. 352 """ 353 if np.any(indices_orig != indices_clipped): 354 raise IERSRangeError('(some) times are outside of range covered ' 355 'by IERS table.') 356 357 def _interpolate(self, jd1, jd2, columns, source=None): 358 mjd, utc = self.mjd_utc(jd1, jd2) 359 # enforce array 360 is_scalar = not hasattr(mjd, '__array__') or mjd.ndim == 0 361 if is_scalar: 362 mjd = np.array([mjd]) 363 utc = np.array([utc]) 364 elif mjd.size == 0: 365 # Short-cut empty input. 366 return np.array([]) 367 368 self._refresh_table_as_needed(mjd) 369 370 # For typical format, will always find a match (since MJD are integer) 371 # hence, important to define which side we will be; this ensures 372 # self['MJD'][i-1]<=mjd<self['MJD'][i] 373 i = np.searchsorted(self['MJD'].value, mjd, side='right') 374 375 # Get index to MJD at or just below given mjd, clipping to ensure we 376 # stay in range of table (status will be set below for those outside) 377 i1 = np.clip(i, 1, len(self) - 1) 378 i0 = i1 - 1 379 mjd_0, mjd_1 = self['MJD'][i0].value, self['MJD'][i1].value 380 results = [] 381 for column in columns: 382 val_0, val_1 = self[column][i0], self[column][i1] 383 d_val = val_1 - val_0 384 if column == 'UT1_UTC': 385 # Check & correct for possible leap second (correcting diff., 386 # not 1st point, since jump can only happen right at 2nd point) 387 d_val -= d_val.round() 388 # Linearly interpolate (which is what TEMPO does for UT1-UTC, but 389 # may want to follow IERS gazette #13 for more precise 390 # interpolation and correction for tidal effects; 391 # https://maia.usno.navy.mil/iers-gaz13) 392 val = val_0 + (mjd - mjd_0 + utc) / (mjd_1 - mjd_0) * d_val 393 394 # Do not extrapolate outside range, instead just propagate last values. 395 val[i == 0] = self[column][0] 396 val[i == len(self)] = self[column][-1] 397 398 if is_scalar: 399 val = val[0] 400 401 results.append(val) 402 403 if source: 404 # Set status to source, using the routine passed in. 405 status = source(i1) 406 # Check for out of range 407 status[i == 0] = TIME_BEFORE_IERS_RANGE 408 status[i == len(self)] = TIME_BEYOND_IERS_RANGE 409 if is_scalar: 410 status = status[0] 411 results.append(status) 412 return results 413 else: 414 self._check_interpolate_indices(i1, i, np.max(mjd)) 415 return results[0] if len(results) == 1 else results 416 417 def _refresh_table_as_needed(self, mjd): 418 """ 419 Potentially update the IERS table in place depending on the requested 420 time values in ``mdj`` and the time span of the table. The base behavior 421 is not to update the table. ``IERS_Auto`` overrides this method. 422 """ 423 pass 424 425 def ut1_utc_source(self, i): 426 """Source for UT1-UTC. To be overridden by subclass.""" 427 return np.zeros_like(i) 428 429 def dcip_source(self, i): 430 """Source for CIP correction. To be overridden by subclass.""" 431 return np.zeros_like(i) 432 433 def pm_source(self, i): 434 """Source for polar motion. To be overridden by subclass.""" 435 return np.zeros_like(i) 436 437 @property 438 def time_now(self): 439 """ 440 Property to provide the current time, but also allow for explicitly setting 441 the _time_now attribute for testing purposes. 442 """ 443 try: 444 return self._time_now 445 except Exception: 446 return Time.now() 447 448 def _convert_col_for_table(self, col): 449 # Fill masked columns with units to avoid dropped-mask warnings 450 # when converting to Quantity. 451 # TODO: Once we support masked quantities, we can drop this and 452 # in the code below replace b_bad with table['UT1_UTC_B'].mask, etc. 453 if (getattr(col, 'unit', None) is not None and 454 isinstance(col, MaskedColumn)): 455 col = col.filled(np.nan) 456 457 return super()._convert_col_for_table(col) 458 459 460class IERS_A(IERS): 461 """IERS Table class targeted to IERS A, provided by USNO. 462 463 These include rapid turnaround and predicted times. 464 See https://datacenter.iers.org/eop.php 465 466 Notes 467 ----- 468 The IERS A file is not part of astropy. It can be downloaded from 469 ``iers.IERS_A_URL`` or ``iers.IERS_A_URL_MIRROR``. See ``iers.__doc__`` 470 for instructions on use in ``Time``, etc. 471 """ 472 473 iers_table = None 474 475 @classmethod 476 def _combine_a_b_columns(cls, iers_a): 477 """ 478 Return a new table with appropriate combination of IERS_A and B columns. 479 """ 480 # IERS A has some rows at the end that hold nothing but dates & MJD 481 # presumably to be filled later. Exclude those a priori -- there 482 # should at least be a predicted UT1-UTC and PM! 483 table = iers_a[np.isfinite(iers_a['UT1_UTC_A']) & 484 (iers_a['PolPMFlag_A'] != '')] 485 486 # This does nothing for IERS_A, but allows IERS_Auto to ensure the 487 # IERS B values in the table are consistent with the true ones. 488 table = cls._substitute_iers_b(table) 489 490 # Combine A and B columns, using B where possible. 491 b_bad = np.isnan(table['UT1_UTC_B']) 492 table['UT1_UTC'] = np.where(b_bad, table['UT1_UTC_A'], table['UT1_UTC_B']) 493 table['UT1Flag'] = np.where(b_bad, table['UT1Flag_A'], 'B') 494 # Repeat for polar motions. 495 b_bad = np.isnan(table['PM_X_B']) | np.isnan(table['PM_Y_B']) 496 table['PM_x'] = np.where(b_bad, table['PM_x_A'], table['PM_X_B']) 497 table['PM_y'] = np.where(b_bad, table['PM_y_A'], table['PM_Y_B']) 498 table['PolPMFlag'] = np.where(b_bad, table['PolPMFlag_A'], 'B') 499 500 b_bad = np.isnan(table['dX_2000A_B']) | np.isnan(table['dY_2000A_B']) 501 table['dX_2000A'] = np.where(b_bad, table['dX_2000A_A'], table['dX_2000A_B']) 502 table['dY_2000A'] = np.where(b_bad, table['dY_2000A_A'], table['dY_2000A_B']) 503 table['NutFlag'] = np.where(b_bad, table['NutFlag_A'], 'B') 504 505 # Get the table index for the first row that has predictive values 506 # PolPMFlag_A IERS (I) or Prediction (P) flag for 507 # Bull. A polar motion values 508 # UT1Flag_A IERS (I) or Prediction (P) flag for 509 # Bull. A UT1-UTC values 510 # Since only 'P' and 'I' are possible and 'P' is guaranteed to come 511 # after 'I', we can use searchsorted for 100 times speed up over 512 # finding the first index where the flag equals 'P'. 513 p_index = min(np.searchsorted(table['UT1Flag_A'], 'P'), 514 np.searchsorted(table['PolPMFlag_A'], 'P')) 515 table.meta['predictive_index'] = p_index 516 table.meta['predictive_mjd'] = table['MJD'][p_index].value 517 518 return table 519 520 @classmethod 521 def _substitute_iers_b(cls, table): 522 # See documentation in IERS_Auto. 523 return table 524 525 @classmethod 526 def read(cls, file=None, readme=None): 527 """Read IERS-A table from a finals2000a.* file provided by USNO. 528 529 Parameters 530 ---------- 531 file : str 532 full path to ascii file holding IERS-A data. 533 Defaults to ``iers.IERS_A_FILE``. 534 readme : str 535 full path to ascii file holding CDS-style readme. 536 Defaults to package version, ``iers.IERS_A_README``. 537 538 Returns 539 ------- 540 ``IERS_A`` class instance 541 """ 542 if file is None: 543 file = IERS_A_FILE 544 if readme is None: 545 readme = IERS_A_README 546 547 iers_a = super().read(file, format='cds', readme=readme) 548 549 # Combine the A and B data for UT1-UTC and PM columns 550 table = cls._combine_a_b_columns(iers_a) 551 table.meta['data_path'] = file 552 table.meta['readme_path'] = readme 553 554 return table 555 556 def ut1_utc_source(self, i): 557 """Set UT1-UTC source flag for entries in IERS table""" 558 ut1flag = self['UT1Flag'][i] 559 source = np.ones_like(i) * FROM_IERS_B 560 source[ut1flag == 'I'] = FROM_IERS_A 561 source[ut1flag == 'P'] = FROM_IERS_A_PREDICTION 562 return source 563 564 def dcip_source(self, i): 565 """Set CIP correction source flag for entries in IERS table""" 566 nutflag = self['NutFlag'][i] 567 source = np.ones_like(i) * FROM_IERS_B 568 source[nutflag == 'I'] = FROM_IERS_A 569 source[nutflag == 'P'] = FROM_IERS_A_PREDICTION 570 return source 571 572 def pm_source(self, i): 573 """Set polar motion source flag for entries in IERS table""" 574 pmflag = self['PolPMFlag'][i] 575 source = np.ones_like(i) * FROM_IERS_B 576 source[pmflag == 'I'] = FROM_IERS_A 577 source[pmflag == 'P'] = FROM_IERS_A_PREDICTION 578 return source 579 580 581class IERS_B(IERS): 582 """IERS Table class targeted to IERS B, provided by IERS itself. 583 584 These are final values; see https://www.iers.org/IERS/EN/Home/home_node.html 585 586 Notes 587 ----- 588 If the package IERS B file (```iers.IERS_B_FILE``) is out of date, a new 589 version can be downloaded from ``iers.IERS_B_URL``. 590 """ 591 592 iers_table = None 593 594 @classmethod 595 def read(cls, file=None, readme=None, data_start=14): 596 """Read IERS-B table from a eopc04_iau2000.* file provided by IERS. 597 598 Parameters 599 ---------- 600 file : str 601 full path to ascii file holding IERS-B data. 602 Defaults to package version, ``iers.IERS_B_FILE``. 603 readme : str 604 full path to ascii file holding CDS-style readme. 605 Defaults to package version, ``iers.IERS_B_README``. 606 data_start : int 607 starting row. Default is 14, appropriate for standard IERS files. 608 609 Returns 610 ------- 611 ``IERS_B`` class instance 612 """ 613 if file is None: 614 file = IERS_B_FILE 615 if readme is None: 616 readme = IERS_B_README 617 618 table = super().read(file, format='cds', readme=readme, 619 data_start=data_start) 620 621 table.meta['data_path'] = file 622 table.meta['readme_path'] = readme 623 return table 624 625 def ut1_utc_source(self, i): 626 """Set UT1-UTC source flag for entries in IERS table""" 627 return np.ones_like(i) * FROM_IERS_B 628 629 def dcip_source(self, i): 630 """Set CIP correction source flag for entries in IERS table""" 631 return np.ones_like(i) * FROM_IERS_B 632 633 def pm_source(self, i): 634 """Set PM source flag for entries in IERS table""" 635 return np.ones_like(i) * FROM_IERS_B 636 637 638class IERS_Auto(IERS_A): 639 """ 640 Provide most-recent IERS data and automatically handle downloading 641 of updated values as necessary. 642 """ 643 iers_table = None 644 645 @classmethod 646 def open(cls): 647 """If the configuration setting ``astropy.utils.iers.conf.auto_download`` 648 is set to True (default), then open a recent version of the IERS-A 649 table with predictions for UT1-UTC and polar motion out to 650 approximately one year from now. If the available version of this file 651 is older than ``astropy.utils.iers.conf.auto_max_age`` days old 652 (or non-existent) then it will be downloaded over the network and cached. 653 654 If the configuration setting ``astropy.utils.iers.conf.auto_download`` 655 is set to False then ``astropy.utils.iers.IERS()`` is returned. This 656 is normally the IERS-B table that is supplied with astropy. 657 658 On the first call in a session, the table will be memoized (in the 659 ``iers_table`` class attribute), and further calls to ``open`` will 660 return this stored table. 661 662 Returns 663 ------- 664 `~astropy.table.QTable` instance 665 With IERS (Earth rotation) data columns 666 667 """ 668 if not conf.auto_download: 669 cls.iers_table = IERS_B.open() 670 return cls.iers_table 671 672 all_urls = (conf.iers_auto_url, conf.iers_auto_url_mirror) 673 674 if cls.iers_table is not None: 675 676 # If the URL has changed, we need to redownload the file, so we 677 # should ignore the internally cached version. 678 679 if cls.iers_table.meta.get('data_url') in all_urls: 680 return cls.iers_table 681 682 try: 683 filename = download_file(all_urls[0], sources=all_urls, cache=True) 684 except Exception as err: 685 # Issue a warning here, perhaps user is offline. An exception 686 # will be raised downstream when actually trying to interpolate 687 # predictive values. 688 warn(AstropyWarning( 689 f'failed to download {" and ".join(all_urls)}, ' 690 f'using local IERS-B: {err}')) 691 cls.iers_table = IERS_B.open() 692 return cls.iers_table 693 694 cls.iers_table = cls.read(file=filename) 695 cls.iers_table.meta['data_url'] = all_urls[0] 696 697 return cls.iers_table 698 699 def _check_interpolate_indices(self, indices_orig, indices_clipped, max_input_mjd): 700 """Check that the indices from interpolation match those after clipping to the 701 valid table range. The IERS_Auto class is exempted as long as it has 702 sufficiently recent available data so the clipped interpolation is 703 always within the confidence bounds of current Earth rotation 704 knowledge. 705 """ 706 predictive_mjd = self.meta['predictive_mjd'] 707 708 # See explanation in _refresh_table_as_needed for these conditions 709 auto_max_age = _none_to_float(conf.auto_max_age) 710 if (max_input_mjd > predictive_mjd and 711 self.time_now.mjd - predictive_mjd > auto_max_age): 712 raise ValueError(INTERPOLATE_ERROR.format(auto_max_age)) 713 714 def _refresh_table_as_needed(self, mjd): 715 """Potentially update the IERS table in place depending on the requested 716 time values in ``mjd`` and the time span of the table. 717 718 For IERS_Auto the behavior is that the table is refreshed from the IERS 719 server if both the following apply: 720 721 - Any of the requested IERS values are predictive. The IERS-A table 722 contains predictive data out for a year after the available 723 definitive values. 724 - The first predictive values are at least ``conf.auto_max_age days`` old. 725 In other words the IERS-A table was created by IERS long enough 726 ago that it can be considered stale for predictions. 727 """ 728 max_input_mjd = np.max(mjd) 729 now_mjd = self.time_now.mjd 730 731 # IERS-A table contains predictive data out for a year after 732 # the available definitive values. 733 fpi = self.meta['predictive_index'] 734 predictive_mjd = self.meta['predictive_mjd'] 735 736 # Update table in place if necessary 737 auto_max_age = _none_to_float(conf.auto_max_age) 738 739 # If auto_max_age is smaller than IERS update time then repeated downloads may 740 # occur without getting updated values (giving a IERSStaleWarning). 741 if auto_max_age < 10: 742 raise ValueError('IERS auto_max_age configuration value must be larger than 10 days') 743 744 if (max_input_mjd > predictive_mjd and 745 (now_mjd - predictive_mjd) > auto_max_age): 746 747 all_urls = (conf.iers_auto_url, conf.iers_auto_url_mirror) 748 749 # Get the latest version 750 try: 751 filename = download_file( 752 all_urls[0], sources=all_urls, cache="update") 753 except Exception as err: 754 # Issue a warning here, perhaps user is offline. An exception 755 # will be raised downstream when actually trying to interpolate 756 # predictive values. 757 warn(AstropyWarning( 758 f'failed to download {" and ".join(all_urls)}: {err}.\n' 759 'A coordinate or time-related ' 760 'calculation might be compromised or fail because the dates are ' 761 'not covered by the available IERS file. See the ' 762 '"IERS data access" section of the astropy documentation ' 763 'for additional information on working offline.')) 764 return 765 766 new_table = self.__class__.read(file=filename) 767 new_table.meta['data_url'] = str(all_urls[0]) 768 769 # New table has new values? 770 if new_table['MJD'][-1] > self['MJD'][-1]: 771 # Replace *replace* current values from the first predictive index through 772 # the end of the current table. This replacement is much faster than just 773 # deleting all rows and then using add_row for the whole duration. 774 new_fpi = np.searchsorted(new_table['MJD'].value, predictive_mjd, side='right') 775 n_replace = len(self) - fpi 776 self[fpi:] = new_table[new_fpi:new_fpi + n_replace] 777 778 # Sanity check for continuity 779 if new_table['MJD'][new_fpi + n_replace] - self['MJD'][-1] != 1.0 * u.d: 780 raise ValueError('unexpected gap in MJD when refreshing IERS table') 781 782 # Now add new rows in place 783 for row in new_table[new_fpi + n_replace:]: 784 self.add_row(row) 785 786 self.meta.update(new_table.meta) 787 else: 788 warn(IERSStaleWarning( 789 'IERS_Auto predictive values are older than {} days but downloading ' 790 'the latest table did not find newer values'.format(conf.auto_max_age))) 791 792 @classmethod 793 def _substitute_iers_b(cls, table): 794 """Substitute IERS B values with those from a real IERS B table. 795 796 IERS-A has IERS-B values included, but for reasons unknown these 797 do not match the latest IERS-B values (see comments in #4436). 798 Here, we use the bundled astropy IERS-B table to overwrite the values 799 in the downloaded IERS-A table. 800 """ 801 iers_b = IERS_B.open() 802 # Substitute IERS-B values for existing B values in IERS-A table 803 mjd_b = table['MJD'][np.isfinite(table['UT1_UTC_B'])] 804 i0 = np.searchsorted(iers_b['MJD'], mjd_b[0], side='left') 805 i1 = np.searchsorted(iers_b['MJD'], mjd_b[-1], side='right') 806 iers_b = iers_b[i0:i1] 807 n_iers_b = len(iers_b) 808 # If there is overlap then replace IERS-A values from available IERS-B 809 if n_iers_b > 0: 810 # Sanity check that we are overwriting the correct values 811 if not u.allclose(table['MJD'][:n_iers_b], iers_b['MJD']): 812 raise ValueError('unexpected mismatch when copying ' 813 'IERS-B values into IERS-A table.') 814 # Finally do the overwrite 815 table['UT1_UTC_B'][:n_iers_b] = iers_b['UT1_UTC'] 816 table['PM_X_B'][:n_iers_b] = iers_b['PM_x'] 817 table['PM_Y_B'][:n_iers_b] = iers_b['PM_y'] 818 table['dX_2000A_B'][:n_iers_b] = iers_b['dX_2000A'] 819 table['dY_2000A_B'][:n_iers_b] = iers_b['dY_2000A'] 820 821 return table 822 823 824class earth_orientation_table(ScienceState): 825 """Default IERS table for Earth rotation and reference systems service. 826 827 These tables are used to calculate the offsets between ``UT1`` and ``UTC`` 828 and for conversion to Earth-based coordinate systems. 829 830 The state itself is an IERS table, as an instance of one of the 831 `~astropy.utils.iers.IERS` classes. The default, the auto-updating 832 `~astropy.utils.iers.IERS_Auto` class, should suffice for most 833 purposes. 834 835 Examples 836 -------- 837 To temporarily use the IERS-B file packaged with astropy:: 838 839 >>> from astropy.utils import iers 840 >>> from astropy.time import Time 841 >>> iers_b = iers.IERS_B.open(iers.IERS_B_FILE) 842 >>> with iers.earth_orientation_table.set(iers_b): 843 ... print(Time('2000-01-01').ut1.isot) 844 2000-01-01T00:00:00.355 845 846 To use the most recent IERS-A file for the whole session:: 847 848 >>> iers_a = iers.IERS_A.open(iers.IERS_A_URL) # doctest: +SKIP 849 >>> iers.earth_orientation_table.set(iers_a) # doctest: +SKIP 850 <ScienceState earth_orientation_table: <IERS_A length=17463>...> 851 852 To go back to the default (of `~astropy.utils.iers.IERS_Auto`):: 853 854 >>> iers.earth_orientation_table.set(None) # doctest: +SKIP 855 <ScienceState earth_orientation_table: <IERS_Auto length=17428>...> 856 """ 857 _value = None 858 859 @classmethod 860 def validate(cls, value): 861 if value is None: 862 value = IERS_Auto.open() 863 if not isinstance(value, IERS): 864 raise ValueError("earth_orientation_table requires an IERS Table.") 865 return value 866 867 868class LeapSeconds(QTable): 869 """Leap seconds class, holding TAI-UTC differences. 870 871 The table should hold columns 'year', 'month', 'tai_utc'. 872 873 Methods are provided to initialize the table from IERS ``Leap_Second.dat``, 874 IETF/ntp ``leap-seconds.list``, or built-in ERFA/SOFA, and to update the 875 list used by ERFA. 876 877 Notes 878 ----- 879 Astropy has a built-in ``iers.IERS_LEAP_SECONDS_FILE``. Up to date versions 880 can be downloaded from ``iers.IERS_LEAP_SECONDS_URL`` or 881 ``iers.LEAP_SECONDS_LIST_URL``. Many systems also store a version 882 of ``leap-seconds.list`` for use with ``ntp`` (e.g., on Debian/Ubuntu 883 systems, ``/usr/share/zoneinfo/leap-seconds.list``). 884 885 To prevent querying internet resources if the available local leap second 886 file(s) are out of date, set ``iers.conf.auto_download = False``. This 887 must be done prior to performing any ``Time`` scale transformations related 888 to UTC (e.g. converting from UTC to TAI). 889 """ 890 # Note: Time instances in this class should use scale='tai' to avoid 891 # needing leap seconds in their creation or interpretation. 892 893 _re_expires = re.compile(r'^#.*File expires on[:\s]+(\d+\s\w+\s\d+)\s*$') 894 _expires = None 895 _auto_open_files = ['erfa', 896 IERS_LEAP_SECOND_FILE, 897 'system_leap_second_file', 898 'iers_leap_second_auto_url', 899 'ietf_leap_second_auto_url'] 900 """Files or conf attributes to try in auto_open.""" 901 902 @classmethod 903 def open(cls, file=None, cache=False): 904 """Open a leap-second list. 905 906 Parameters 907 ---------- 908 file : path-like or None 909 Full local or network path to the file holding leap-second data, 910 for passing on to the various ``from_`` class methods. 911 If 'erfa', return the data used by the ERFA library. 912 If `None`, use default locations from file and configuration to 913 find a table that is not expired. 914 cache : bool 915 Whether to use cache. Defaults to False, since leap-second files 916 are regularly updated. 917 918 Returns 919 ------- 920 leap_seconds : `~astropy.utils.iers.LeapSeconds` 921 Table with 'year', 'month', and 'tai_utc' columns, plus possibly 922 others. 923 924 Notes 925 ----- 926 Bulletin C is released about 10 days after a possible leap second is 927 introduced, i.e., mid-January or mid-July. Expiration days are thus 928 generally at least 150 days after the present. For the auto-loading, 929 a list comprised of the table shipped with astropy, and files and 930 URLs in `~astropy.utils.iers.Conf` are tried, returning the first 931 that is sufficiently new, or the newest among them all. 932 """ 933 if file is None: 934 return cls.auto_open() 935 936 if file.lower() == 'erfa': 937 return cls.from_erfa() 938 939 if urlparse(file).netloc: 940 file = download_file(file, cache=cache) 941 942 # Just try both reading methods. 943 try: 944 return cls.from_iers_leap_seconds(file) 945 except Exception: 946 return cls.from_leap_seconds_list(file) 947 948 @staticmethod 949 def _today(): 950 # Get current day in scale='tai' without going through a scale change 951 # (so we do not need leap seconds). 952 s = '{0.year:04d}-{0.month:02d}-{0.day:02d}'.format(datetime.utcnow()) 953 return Time(s, scale='tai', format='iso', out_subfmt='date') 954 955 @classmethod 956 def auto_open(cls, files=None): 957 """Attempt to get an up-to-date leap-second list. 958 959 The routine will try the files in sequence until it finds one 960 whose expiration date is "good enough" (see below). If none 961 are good enough, it returns the one with the most recent expiration 962 date, warning if that file is expired. 963 964 For remote files that are cached already, the cached file is tried 965 first before attempting to retrieve it again. 966 967 Parameters 968 ---------- 969 files : list of path-like, optional 970 List of files/URLs to attempt to open. By default, uses 971 ``cls._auto_open_files``. 972 973 Returns 974 ------- 975 leap_seconds : `~astropy.utils.iers.LeapSeconds` 976 Up to date leap-second table 977 978 Notes 979 ----- 980 Bulletin C is released about 10 days after a possible leap second is 981 introduced, i.e., mid-January or mid-July. Expiration days are thus 982 generally at least 150 days after the present. We look for a file 983 that expires more than 180 - `~astropy.utils.iers.Conf.auto_max_age` 984 after the present. 985 """ 986 good_enough = cls._today() + TimeDelta(180-_none_to_float(conf.auto_max_age), 987 format='jd') 988 989 if files is None: 990 # Basic files to go over (entries in _auto_open_files can be 991 # configuration items, which we want to be sure are up to date). 992 files = [getattr(conf, f, f) for f in cls._auto_open_files] 993 994 # Remove empty entries. 995 files = [f for f in files if f] 996 997 # Our trials start with normal files and remote ones that are 998 # already in cache. The bools here indicate that the cache 999 # should be used. 1000 trials = [(f, True) for f in files 1001 if not urlparse(f).netloc or is_url_in_cache(f)] 1002 # If we are allowed to download, we try downloading new versions 1003 # if none of the above worked. 1004 if conf.auto_download: 1005 trials += [(f, False) for f in files if urlparse(f).netloc] 1006 1007 self = None 1008 err_list = [] 1009 # Go through all entries, and return the first one that 1010 # is not expired, or the most up to date one. 1011 for f, allow_cache in trials: 1012 if not allow_cache: 1013 clear_download_cache(f) 1014 1015 try: 1016 trial = cls.open(f, cache=True) 1017 except Exception as exc: 1018 err_list.append(exc) 1019 continue 1020 1021 if self is None or trial.expires > self.expires: 1022 self = trial 1023 self.meta['data_url'] = str(f) 1024 if self.expires > good_enough: 1025 break 1026 1027 if self is None: 1028 raise ValueError('none of the files could be read. The ' 1029 'following errors were raised:\n' + str(err_list)) 1030 1031 if self.expires < self._today(): 1032 warn('leap-second file is expired.', IERSStaleWarning) 1033 1034 return self 1035 1036 @property 1037 def expires(self): 1038 """The limit of validity of the table.""" 1039 return self._expires 1040 1041 @classmethod 1042 def _read_leap_seconds(cls, file, **kwargs): 1043 """Read a file, identifying expiration by matching 'File expires'""" 1044 expires = None 1045 # Find expiration date. 1046 with get_readable_fileobj(file) as fh: 1047 lines = fh.readlines() 1048 for line in lines: 1049 match = cls._re_expires.match(line) 1050 if match: 1051 day, month, year = match.groups()[0].split() 1052 month_nb = MONTH_ABBR.index(month[:3]) + 1 1053 expires = Time(f'{year}-{month_nb:02d}-{day}', 1054 scale='tai', out_subfmt='date') 1055 break 1056 else: 1057 raise ValueError(f'did not find expiration date in {file}') 1058 1059 self = cls.read(lines, format='ascii.no_header', **kwargs) 1060 self._expires = expires 1061 return self 1062 1063 @classmethod 1064 def from_iers_leap_seconds(cls, file=IERS_LEAP_SECOND_FILE): 1065 """Create a table from a file like the IERS ``Leap_Second.dat``. 1066 1067 Parameters 1068 ---------- 1069 file : path-like, optional 1070 Full local or network path to the file holding leap-second data 1071 in a format consistent with that used by IERS. By default, uses 1072 ``iers.IERS_LEAP_SECOND_FILE``. 1073 1074 Notes 1075 ----- 1076 The file *must* contain the expiration date in a comment line, like 1077 '# File expires on 28 June 2020' 1078 """ 1079 return cls._read_leap_seconds( 1080 file, names=['mjd', 'day', 'month', 'year', 'tai_utc']) 1081 1082 @classmethod 1083 def from_leap_seconds_list(cls, file): 1084 """Create a table from a file like the IETF ``leap-seconds.list``. 1085 1086 Parameters 1087 ---------- 1088 file : path-like, optional 1089 Full local or network path to the file holding leap-second data 1090 in a format consistent with that used by IETF. Up to date versions 1091 can be retrieved from ``iers.IETF_LEAP_SECOND_URL``. 1092 1093 Notes 1094 ----- 1095 The file *must* contain the expiration date in a comment line, like 1096 '# File expires on: 28 June 2020' 1097 """ 1098 from astropy.io.ascii import convert_numpy # Here to avoid circular import 1099 1100 names = ['ntp_seconds', 'tai_utc', 'comment', 'day', 'month', 'year'] 1101 # Note: ntp_seconds does not fit in 32 bit, so causes problems on 1102 # 32-bit systems without the np.int64 converter. 1103 self = cls._read_leap_seconds( 1104 file, names=names, include_names=names[:2], 1105 converters={'ntp_seconds': [convert_numpy(np.int64)]}) 1106 self['mjd'] = (self['ntp_seconds']/86400 + 15020).round() 1107 # Note: cannot use Time.ymdhms, since that might require leap seconds. 1108 isot = Time(self['mjd'], format='mjd', scale='tai').isot 1109 ymd = np.array([[int(part) for part in t.partition('T')[0].split('-')] 1110 for t in isot]) 1111 self['year'], self['month'], self['day'] = ymd.T 1112 return self 1113 1114 @classmethod 1115 def from_erfa(cls, built_in=False): 1116 """Create table from the leap-second list in ERFA. 1117 1118 Parameters 1119 ---------- 1120 built_in : bool 1121 If `False` (default), retrieve the list currently used by ERFA, 1122 which may have been updated. If `True`, retrieve the list shipped 1123 with erfa. 1124 """ 1125 current = cls(erfa.leap_seconds.get()) 1126 current._expires = Time('{0.year:04d}-{0.month:02d}-{0.day:02d}' 1127 .format(erfa.leap_seconds.expires), 1128 scale='tai') 1129 if not built_in: 1130 return current 1131 1132 try: 1133 erfa.leap_seconds.set(None) # reset to defaults 1134 return cls.from_erfa(built_in=False) 1135 finally: 1136 erfa.leap_seconds.set(current) 1137 1138 def update_erfa_leap_seconds(self, initialize_erfa=False): 1139 """Add any leap seconds not already present to the ERFA table. 1140 1141 This method matches leap seconds with those present in the ERFA table, 1142 and extends the latter as necessary. 1143 1144 Parameters 1145 ---------- 1146 initialize_erfa : bool, or 'only', or 'empty' 1147 Initialize the ERFA leap second table to its built-in value before 1148 trying to expand it. This is generally not needed but can help 1149 in case it somehow got corrupted. If equal to 'only', the ERFA 1150 table is reinitialized and no attempt it made to update it. 1151 If 'empty', the leap second table is emptied before updating, i.e., 1152 it is overwritten altogether (note that this may break things in 1153 surprising ways, as most leap second tables do not include pre-1970 1154 pseudo leap-seconds; you were warned). 1155 1156 Returns 1157 ------- 1158 n_update : int 1159 Number of items updated. 1160 1161 Raises 1162 ------ 1163 ValueError 1164 If the leap seconds in the table are not on 1st of January or July, 1165 or if the matches are inconsistent. This would normally suggest 1166 a corrupted leap second table, but might also indicate that the 1167 ERFA table was corrupted. If needed, the ERFA table can be reset 1168 by calling this method with an appropriate value for 1169 ``initialize_erfa``. 1170 """ 1171 if initialize_erfa == 'empty': 1172 # Initialize to empty and update is the same as overwrite. 1173 erfa.leap_seconds.set(self) 1174 return len(self) 1175 1176 if initialize_erfa: 1177 erfa.leap_seconds.set() 1178 if initialize_erfa == 'only': 1179 return 0 1180 1181 return erfa.leap_seconds.update(self) 1182