casa6
Public
Commits
| 1 - | # main module for retreiving JPL-Horizons ephemeris data and |
| 2 - | # convert it to CASA readable format (i.e. CASA table ) |
| 3 - | import numpy as np |
| 1 + | ''' |
| 2 + | main module for retreiving JPL-Horizons ephemeris data and |
| 3 + | convert it to CASA readable format (i.e. CASA table ) |
| 4 + | ''' |
| 4 5 | import os |
| 5 6 | import time |
| 6 7 | from math import sqrt, sin, cos, log |
| 8 + | import numpy as np |
| 7 9 | |
| 8 10 | from casatools import table, quanta, measures |
| 9 11 | from casatasks import casalog |
| 10 12 | |
| 11 13 | _tb = table() |
| 12 14 | _qa = quanta() |
| 13 15 | _me = measures() |
| 14 16 | |
| 15 17 | # debug |
| 16 18 | _debug = False |
| 26 28 | format |
| 27 29 | stoptime: stop (end) time of the ephemeris data in the same format as |
| 28 30 | starttime |
| 29 31 | incr: time increment (interval) of the ephemeris data |
| 30 32 | outtable: output CASA table name |
| 31 33 | asis: a toggle for additinal check to resolve object name |
| 32 34 | rawdatafile: raw query result file name (optional) |
| 33 35 | """ |
| 34 36 | |
| 35 37 | # commented ones are not currently supported in setjy |
| 36 - | # Juno's table is exist in the data repo but not supported in setjy |
| 38 + | # Juno's table is exist in the data repo but not supported in setjy |
| 37 39 | asteroids = {'ceres': '1;', |
| 38 40 | 'pallas': '2;', |
| 39 41 | 'juno': '3;', # large crater and temperature varies |
| 40 42 | 'vesta': '4;', |
| 41 43 | 'lutetia': '21;' |
| 42 44 | } |
| 43 45 | planets_and_moons = {'sun': '10', |
| 44 46 | 'mercury': '199', |
| 45 47 | 'venus': '299', |
| 46 48 | 'moon': '301', |
| 65 67 | # it will skip the check. |
| 66 68 | # For small bodies such as comets, objectname must include 'DES=' following designated |
| 67 69 | # object name or id (SPK ID). Based on the JPL-Horizons documentaiton, for IAU ID, 'DES=' |
| 68 70 | # can be omitted. https://ssd.jpl.nasa.gov/horizons/app.html#command |
| 69 71 | start_time = None |
| 70 72 | stop_time = None |
| 71 73 | step_size = None |
| 72 74 | if not asis: |
| 73 75 | if not objectname.lower() in known_objects: |
| 74 76 | raise ValueError( |
| 75 - | "%s is not in the known object list for CASA. To skip this check set asis=True" % objectname) |
| 77 + | f"{objectname} is not in the known object list for CASA. To skip this check set asis=True") |
| 76 78 | else: |
| 77 79 | target = known_objects[objectname.lower()] |
| 78 80 | else: |
| 79 81 | target = objectname |
| 80 82 | try: |
| 81 83 | if starttime.startswith('JD') or starttime.startswith('MJD'): |
| 82 84 | start_time = starttime |
| 83 85 | stop_time = stoptime |
| 84 86 | else: |
| 85 87 | start_time = _qa.time(starttime, form='ymd') |
| 125 127 | from urllib.error import URLError |
| 126 128 | import socket |
| 127 129 | |
| 128 130 | context = ssl.create_default_context(cafile=certifi.where()) |
| 129 131 | urlbase = 'https://ssd.jpl.nasa.gov/api/horizons.api' |
| 130 132 | data = None |
| 131 133 | |
| 132 134 | values = {'format': 'json', |
| 133 135 | 'EPHEM_TYPE': 'OBSERVER', |
| 134 136 | 'OBJ_DATA': 'YES', |
| 135 - | 'COMMAND': "'{}'".format(target), |
| 137 + | 'COMMAND': f"'{target}'", |
| 136 138 | 'START_TIME': starttime, |
| 137 139 | 'STOP_TIME': stoptime, |
| 138 140 | 'STEP_SIZE': stepsize, |
| 139 141 | 'CENTER': '500@399', |
| 140 - | 'QUANTITIES': "'{}'".format(quantities), |
| 142 + | 'QUANTITIES': f"'{quantities}'", |
| 141 143 | 'ANG_FORMAT': ang_format} |
| 142 144 | |
| 143 145 | pardata = urlencode(values, doseq=True, encoding='utf-8') |
| 144 146 | params = pardata.encode('ascii') |
| 145 147 | #print("params=", params) |
| 146 148 | |
| 147 149 | req = Request(urlbase, params) |
| 148 150 | datastr = None |
| 149 151 | try: |
| 150 152 | with urlopen(req, context=context, timeout=10) as response: |
| 151 153 | datastr = response.read().decode() |
| 152 154 | except URLError as e: |
| 153 155 | if hasattr(e, 'reason'): |
| 154 - | msg = 'URLError: Failed to reach URL={0} (reason: {1})'.format(urlbase, e.reason) |
| 156 + | msg = f'URLError: Failed to reach URL={urlbase} (reason: {e.reason})' |
| 155 157 | casalog.post(msg, 'WARN') |
| 156 158 | if hasattr(e, 'code'): |
| 157 - | msg = 'URLError: Couldn\'t fullfill the request to {0} (code: {1})'.format(urlbase, e.code) |
| 159 + | msg = f'URLError: Couldn\'t fullfill the request to {urlbase} (code: {e.code})' |
| 158 160 | casalog.post(msg, 'WARN') |
| 159 161 | except socket.timeout as e: |
| 160 - | msg = 'Failed to reach URL={0}. Socket timeout {1}'.format(urlbase, e) |
| 162 + | msg = f'Failed to reach URL={urlbase}. Socket timeout {e}' |
| 161 163 | casalog.post(msg, 'WARN') |
| 162 164 | |
| 163 165 | status = response.getcode() |
| 164 166 | if datastr: |
| 165 167 | try: |
| 166 168 | data = ast.literal_eval(datastr) |
| 167 169 | except RuntimeError as e: |
| 168 170 | casalog.post(e, 'SEVERE') |
| 169 171 | if status == 200: |
| 170 172 | # |
| 171 173 | # If the ephemeris data file was generated, write it to the output file: |
| 172 174 | if rawdatafile != "": |
| 173 175 | if os.path.exists(rawdatafile): |
| 174 176 | os.remove(rawdatafile) |
| 175 177 | if "result" in data: |
| 176 178 | with open(rawdatafile, "w") as f: |
| 177 179 | f.write(data["result"]) |
| 178 180 | # Otherwise, the ephemeris file was not generated so output an error |
| 179 181 | else: |
| 180 182 | casalog.post("ERROR: No data found. Ephemeris data file not generated", 'WARN') |
| 181 183 | else: |
| 182 - | raise RuntimeError('Could not retrieve the data. Error code:{}:{}'.format(status, response.msg)) |
| 184 + | raise RuntimeError(f'Could not retrieve the data. Error code:{status}:{response.msg}') |
| 183 185 | else: |
| 184 186 | data = None |
| 185 187 | return data |
| 186 188 | |
| 187 189 | |
| 188 190 | def tocasatb(indata, outtable): |
| 189 191 | """ |
| 190 192 | convert a JPL-Horizons query results to a CASA table |
| 191 193 | indata: either query data ('results') or file name that contains the result data |
| 192 194 | outtable: output CASA table name |
| 243 245 | 'unit': 'deg'} |
| 244 246 | } |
| 245 247 | |
| 246 248 | # do in a two-step |
| 247 249 | # read the original query result dictionary containing ephemeris data and save the data part |
| 248 250 | # to a temporary text file. Then re-read the temporary file to a casa table |
| 249 251 | # with the approriate data format that setjy and measure expect. |
| 250 252 | tempfname = 'temp_ephem_'+str(os.getpid())+'.dat' |
| 251 253 | tempconvfname = 'temp_ephem_conv_'+str(os.getpid())+'.dat' |
| 252 254 | try: |
| 255 + | exeedthelinelimit=None |
| 253 256 | # Scan the original data |
| 254 - | if type(indata) == dict and 'result' in indata: |
| 257 + | if isinstance(indata, dict) and 'result' in indata: |
| 255 258 | #print("ephem data dict") |
| 256 259 | ephemdata = indata['result'] |
| 257 - | elif type(indata) == str: |
| 260 + | elif isinstance(indata, str): |
| 258 261 | if os.path.exists(indata): |
| 259 262 | with open(indata, 'r') as infile: |
| 260 263 | ephemdata = infile.readlines() |
| 261 264 | # the relevant information in the main header section is extracted as |
| 262 265 | # it is read line by line. The ephemeris data is stored in the separate text file |
| 263 266 | # to be further processed in subsequent steps. |
| 264 267 | with open(tempfname, 'w') as outfile: |
| 265 268 | # Some initialization |
| 266 269 | headerdict = {} |
| 267 270 | # jplhorizonsdataIdFound = False |
| 268 271 | datalines = 0 |
| 269 272 | readthedata = False |
| 270 273 | startmjd = None |
| 274 + | endmjd = None |
| 271 275 | incolnames = None |
| 272 276 | # multiple entries (in different units) may exit for orb. per. |
| 273 277 | foundorbper = False |
| 274 278 | ### |
| 275 279 | lcnt = 0 |
| 276 280 | # for lnum, line in enumerate(infile): |
| 277 281 | for lnum, line in enumerate(ephemdata.split('\n')): |
| 278 282 | # JPL-Horizons data should contain this line at the beginning |
| 279 283 | if re.search(r'JPL/HORIZONS', line): |
| 280 284 | # jplhorizondataIdFound = True |
| 304 308 | m = re.match(r'^[>\s]*Target body name:\s+\d*\s*(\w+)', line) |
| 305 309 | if m: |
| 306 310 | headerdict['NAME'] = m[1] |
| 307 311 | # start time (of the requested time range) |
| 308 312 | elif re.search(r'Start time', line): |
| 309 313 | m = re.match(r'^[>\s]*Start time\s+\S+\s+\S+\s+(\S+)\s+(\S+)\s+(\w+)', line) |
| 310 314 | if m: |
| 311 315 | startmjd = _qa.totime(m[1] + '/' + m[2]) |
| 312 316 | #--This info will not be used but left here since it might be useful fo debugging. |
| 313 317 | # # end time (of the requested time range) |
| 314 - | # elif re.search(r'End time', line): |
| 315 - | # m = re.match(r'^[>\s]*End time\s+\S+\s+\S+\s+(\S+)\s+(\S+)\s+(\w+)', line) |
| 316 - | # if m: |
| 317 - | # endmjd = _qa.totime(m[1] + '/' + m[2]) |
| 318 + | elif re.search(r'Stop time', line): |
| 319 + | m = re.match(r'^[>\s]*Stop time\s+\S+\s+\S+\s+(\S+)\s+(\S+)\s+(\w+)', line) |
| 320 + | if m: |
| 321 + | endmjd = _qa.totime(m[1] + '/' + m[2]) |
| 318 322 | # date increment |
| 319 323 | elif re.search(r'Step-size', line): |
| 320 324 | m = re.match(r'^[>\s]*Step-size\s+\S+\s+(\S+)\s+(\w+)', line) |
| 321 325 | if m: |
| 322 326 | unit = m[2] |
| 323 327 | if unit == 'minutes': |
| 324 328 | # this is the default unit returned by the JPL-Horizons |
| 325 329 | theunit = 'min' |
| 326 330 | elif unit == 'hours': |
| 327 331 | theunit = 'h' |
| 328 332 | elif unit == 'days': |
| 329 333 | theunit = 'd' |
| 334 + | elif unit == 'steps': |
| 335 + | theunit = 'steps' |
| 336 + | elif unit == 'calendar': |
| 337 + | raise RuntimeError('Unit of Step-size in calendar month or year is not supported.') |
| 330 338 | else: |
| 331 - | raise RuntimeError('Unit of Step-size, %s is unrecognized' % unit) |
| 339 + | raise RuntimeError(f'Unit of Step-size, {unit} is unrecognized') |
| 332 340 | if theunit == 'd': |
| 333 341 | dmjd = m[1] |
| 342 + | elif theunit == 'steps': |
| 343 + | # print('endtime=',endmjd['value']) |
| 344 + | dmjd = (endmjd['value'] - startmjd['value'])/int(m[1]) |
| 334 345 | else: |
| 335 - | dmjd = _qa.convert(_qa.totime(m[1] + theunit), 'd') |
| 336 - | headerdict['dMJD'] = dmjd['value'] |
| 346 + | dmjd = _qa.convert(_qa.totime(m[1] + theunit), 'd')['value'] |
| 347 + | headerdict['dMJD'] = dmjd |
| 337 348 | if startmjd is not None: # start mjd should be available before step-size line |
| 338 349 | # MJD0 = firstMJD - dMJD (as defined casacore MeasComet documentation) |
| 339 - | headerdict['MJD0'] = startmjd['value'] - dmjd['value'] |
| 350 + | headerdict['MJD0'] = startmjd['value'] - dmjd |
| 340 351 | elif re.search(r'Center geodetic', line): |
| 341 352 | m = re.match(r'^[>\s]*Center geodetic\s*: ([-+0-9.]+,\s*[-+0-9.]+,\s*[-+0-9.]+)', line) |
| 342 353 | if m: |
| 343 354 | long_lat_alt = m[1].split(',') |
| 344 355 | headerdict['GeoLong'] = float(long_lat_alt[0]) |
| 345 356 | headerdict['GeoLat'] = float(long_lat_alt[1]) |
| 346 357 | headerdict['GeoDist'] = float(long_lat_alt[2]) |
| 347 358 | # obs location |
| 348 359 | elif re.search(r'Center-site name', line): |
| 349 360 | m = re.match(r'^[>\s]*Center-site name:\s+(\w+)', line) |
| 357 368 | matchlist = m.groups() |
| 358 369 | radiiarr = np.zeros(3) |
| 359 370 | if len(matchlist)==2: |
| 360 371 | if m[2] is None: |
| 361 372 | radiiarr = np.asarray(np.array(m[1].split(' x ')), dtype=np.float64) |
| 362 373 | headerdict['radii'] = {'unit': 'km', 'value': radiiarr} |
| 363 374 | elif m[1] is None: |
| 364 375 | radiiarr = np.array([m[2],m[2],m[2]], dtype=np.float64) |
| 365 376 | headerdict['radii'] = {'unit': 'km', 'value': radiiarr} |
| 366 377 | else: |
| 367 - | casaloog.post("Unexpected number or matches for Target radii:{} (expected 2)".format(m.groups), 'WARN') |
| 378 + | casalog.post(f"Unexpected number or matches for Target radii:{m.groups} (expected 2)", 'WARN') |
| 368 379 | #rotational period (few pattens seem to exist) |
| 369 380 | elif re.search(r'rot. period|Rotational period', line): |
| 370 381 | m = re.search(r'rot. period\s+\S*=\s*([0-9.]+h\s*[0-9.]+m\s*[0-9.]+\s*s)|' |
| 371 382 | 'rot. period\s+\S*=\s*([0-9.]+)(?:\s*\+-[0-9.]+)?\s*([dh])|' |
| 372 383 | 'Rotational period\s*=\s+Synchronous', line) |
| 373 384 | if m: |
| 374 385 | if m[0].find('Synchronous') != -1: |
| 375 386 | headerdict['rot_per'] = 'Synchronous' |
| 376 387 | else: |
| 377 388 | if len(m.groups()) == 3: |
| 410 421 | elif re.search(r'Mean [Tt]emperature', line): |
| 411 422 | m = re.search(r'Mean [Tt]emperature\s+\(K\)\s+=\s+([0-9.]+)\s+',line) |
| 412 423 | if m: |
| 413 424 | headerdict['T_mean'] = {'unit':'K', 'value':float(m[1])} |
| 414 425 | # start reading data |
| 415 426 | elif re.search(r'\s*Date__\(UT\)', line): |
| 416 427 | incolnames = line.split() |
| 417 428 | elif re.search(r'\$\$SOE', line): |
| 418 429 | readthedata = True |
| 419 430 | elif re.search(r'\$\$EOE', line): |
| 420 - | endofdata = True |
| 421 431 | readthedata = False |
| 422 432 | elif readthedata: |
| 423 433 | datalines += 1 |
| 424 434 | outfile.write(line + '\n') |
| 435 + | elif re.search(r'Projected output length', line): |
| 436 + | # this message occurs when requested data is too large |
| 437 + | exeedthelinelimit = line |
| 425 438 | else: |
| 426 439 | pass |
| 427 440 | lcnt += 1 |
| 428 441 | if 'radii' in headerdict: |
| 429 442 | radiival = headerdict['radii']['value'] |
| 430 443 | meanrad = _mean_radius(radiival[0], radiival[1], radiival[2]) |
| 431 444 | headerdict['meanrad'] = {'unit': 'km', 'value': meanrad} |
| 432 - | casalog.post(f"Number of data lines={datalines}") |
| 445 + | if datalines == 0: |
| 446 + | casalog.post("No ephemeris data was found", "WARN") |
| 447 + | if exeedthelinelimit is not None: |
| 448 + | raise RuntimeError("Error occur at the query:"+exeedthelinelimit) |
| 449 + | else: |
| 450 + | casalog.post(f"Number of data lines={datalines}") |
| 433 451 | casalog.post(f"Number of all lines in the file={lcnt}") |
| 434 452 | #print("headerdict=", headerdict) |
| 435 453 | # output to a casa table |
| 436 454 | |
| 437 455 | # check the input data columns and stored the order as indices |
| 438 456 | foundncols = 0 |
| 439 457 | indexoffset = 0 |
| 440 458 | colkeys = {} |
| 441 459 | if incolnames is not None: |
| 442 460 | for outcolname in cols: |
| 466 484 | # Dec data col is next to RA data col |
| 467 485 | cols[outcolname]['index'] = cols['RA']['index'] + 1 |
| 468 486 | # add additional offset for index (4 data columns at this point) |
| 469 487 | indexoffset +=1 |
| 470 488 | foundncols += 1 |
| 471 489 | else: |
| 472 490 | if inheadername in incolnames: |
| 473 491 | cols[outcolname]['index'] = incolnames.index(inheadername) + indexoffset |
| 474 492 | foundncols += 1 |
| 475 493 | else: |
| 476 - | casalog.post(f"Cannot find {ihheadername}", "WARN") |
| 494 + | casalog.post(f"Cannot find {inheadername}", "WARN") |
| 477 495 | |
| 478 496 | #print(cols) |
| 479 497 | casalog.post(f"expected n cols = {len(cols)} ") |
| 480 498 | casalog.post(f"foundncols = {foundncols}") |
| 481 499 | if foundncols == len(cols): |
| 482 500 | # Format the data to comply with measure/setjy |
| 483 501 | casalog.post("Found all the required columns") |
| 484 502 | with open(tempconvfname, 'w') as outf, open(tempfname, 'r') as inf: |
| 485 503 | ndata = 0 |
| 486 504 | earliestmjd = None |
| 535 553 | |
| 536 554 | else: |
| 537 555 | missingcols = len(cols) - foundncols |
| 538 556 | casalog.post(r"Missing {missingcols}") |
| 539 557 | |
| 540 558 | # final step: convert to a CASA table |
| 541 559 | dtypes = np.array(['D' for _ in range(len(cols))]) |
| 542 560 | _tb.fromascii(outtable, tempconvfname, sep=' ', columnnames=list(cols.keys()), |
| 543 561 | datatypes=dtypes.tolist()) |
| 544 562 | _tb.done() |
| 563 + | |
| 545 564 | # fill keyword values in the ephem table |
| 546 565 | if os.path.exists(outtable): |
| 547 566 | _fill_keywords_from_dict(headerdict, colkeys, outtable) |
| 548 567 | casalog.post(f"Output is written to a CASA table, {outtable}") |
| 549 568 | else: |
| 550 569 | raise RuntimeError("Error occured. The output table, " + outtable + "is not generated") |
| 551 - | except RuntimeError: |
| 552 - | raise Exception("Error occurred") |
| 570 + | else: #incolnames is None |
| 571 + | raise RuntimeError("No data header was found.") |
| 572 + | |
| 573 + | except RuntimeError as e: |
| 574 + | casalog.post(str(e),"SEVERE") |
| 553 575 | finally: |
| 554 576 | tempfiles = [tempfname, tempconvfname] |
| 555 577 | _clean_up(tempfiles) |
| 556 578 | |
| 557 579 | # This needs to change to read the ephem data from the generated casa table. |
| 558 580 | # Will be called in fill_keywords_from_dict() |
| 559 581 | |
| 560 582 | # Copied from JPLephem_reader2.py for mean_radius calc (when no NP-ang, NP_dist) |
| 561 583 | # Note: this is not fully verified for the correctness but has been used |
| 562 584 | # through JPLephem_reader2 for processing the ephemeris data to a CASA table. |
| 653 675 | _tb.flush() |
| 654 676 | _tb.done() |
| 655 677 | except RuntimeError: |
| 656 678 | casalog.post('Internal error: Cannot add the data in keywords', "WARN") |
| 657 679 | |
| 658 680 | def _clean_up(filelist): |
| 659 681 | """ |
| 660 682 | Clean up the temporary files |
| 661 683 | """ |
| 662 684 | for f in filelist: |
| 663 - | if os.path.exists(f): |
| 664 - | os.remove(f) |
| 685 + | if os.path.exists(f): |
| 686 + | os.remove(f) |