Skip to content

Commit

Permalink
Cleaned up scripts in data/.
Browse files Browse the repository at this point in the history
Cleaned up aop.py, map.py, summary_test_slalib.py
  aop.py, map.py :Use Numpy array operations when possible.
  In above and summary_test_slalib.py:
      Add ".." to path only if import fails.
Cleaned up data/summary_test_slalib.py.
    summary_test_slalib.py: print std, instead of variance.
    data/read_data.get_sla(): return Numpy array.
    data/slalib_convert_hip.py: remove trailing space from format
      string for ecleq, eqecl, eqgal and galeq.
    pytpm/tests/data/slalib_hip_<above>.txt were edited to
      remove the trailing space.
  • Loading branch information
phn committed Aug 25, 2011
1 parent 3af7ca8 commit 4b4c4a7
Show file tree
Hide file tree
Showing 9 changed files with 4,864 additions and 4,882 deletions.
34 changes: 18 additions & 16 deletions data/aop.py
Original file line number Diff line number Diff line change
@@ -1,21 +1,22 @@
from __future__ import print_function
from __future__ import absolute_import
import sys
import math
from scipy import stats
import numpy as np
import os
from read_data import testdatadir
from read_data import get_hipdata, get_sla
from read_data import get_hipdata
from read_data import cat2array

# I want to run these without having to install PyTPM.
sys.path.append("..")
from pytpm import tpm, convert
try:
from pytpm import tpm, convert
except ImportError:
sys.path.append("..")
from pytpm import tpm, convert

# Hipparcos data
# Load HIPPARCOS data and output form running SLALIB aop.
hip_tab = get_hipdata()

# Result from running SLALIB MAP + AOP functions.
tab = np.loadtxt(os.path.join(testdatadir, "slalib_hip_aop.txt"))

# Normalize "alpha" angles to 0 - 360.0
Expand All @@ -36,16 +37,15 @@
utc = tpm.gcal2j(2010, 1, 1) - 0.5 # midnight
tt = tpm.utc2tdb(utc)

# Convert J2000 RA, DEC to Az, EL and ZD at UTC.
# Convert J2000 RA, DEC to Az, EL and ZD at given UTC.
v6o = convert.proper_motion(v6l, tt, tpm.J2000)
v619 = convert.convertv6(v6o, s1=6, s2=19, utc=utc)
cat19 = convert.v62cat(v619, tpm.CJ)
cat19 = cat2array(cat19)

az = [math.degrees(j['alpha']) for j in cat19]
el = [math.degrees(j['delta']) for j in cat19]
az = np.array(az)
zd = 90 - np.array(el)

az = np.degrees(cat19['alpha'])
el = np.degrees(cat19['delta'])
zd = 90.0 - el
# Keep only those objects with ZD < 75.0 degrees.
indx = np.where(zd < 75.0)

Expand All @@ -56,9 +56,10 @@
# Az, El to HA and Dec.
v620 = convert.convertv6(v619, s1=19, s2=20, utc=utc)
cat20 = convert.v62cat(v620, tpm.CJ)
cat20 = cat2array(cat20)

ha = np.array([math.degrees(j['alpha']) for j in cat20])
dec = np.array([math.degrees(j['delta']) for j in cat20])
ha = np.degrees(cat20['alpha'])
dec = np.degrees(cat20['delta'])

# Difference in HA and Dec, using TPM and SLALIB.
ha_diff = np.abs(ha[indx] - ha_sla[indx]) * 3600.0
Expand All @@ -75,10 +76,11 @@
tpm.tpm_data(tstate, tpm.TPM_ALL)
last = tpm.r2d(tstate.last)
ra = last - ha
# Have to normalize to 0 - 360.0
# Have to normalize to 0 - 360.0.
ra = np.array([i if i > 0 else i + 360.0 for i in ra])
ra_diff = np.abs(ra[indx] - ra_sla[indx]) * 3600.0

print("Comparison with SLALIB aop using HIPPARCOS data.")
fs = "{0} {1}\n" + \
"Min: {2:.4f} Max: {3:.4f} \nMean: {4:.4f} Std: {5:.4f}\n"
x = stats.describe(az_diff)
Expand Down
29 changes: 16 additions & 13 deletions data/map.py
Original file line number Diff line number Diff line change
@@ -1,43 +1,46 @@
from __future__ import print_function
from __future__ import absolute_import
import sys
import math
from scipy import stats
import numpy as np
import os
from read_data import testdatadir
from read_data import get_hipdata, get_sla
from read_data import get_hipdata
from read_data import cat2array

# I want to run these without having to install PyTPM.
sys.path.append("..")
from pytpm import tpm, convert


#tab = get_sla("slalib_hip_map.txt")
tab = np.loadtxt(os.path.join(testdatadir, "slalib_hip_map.txt"))
# Load HIPPAROCS data and output from running SLALIB map.
hip_tab = get_hipdata()

tab = np.loadtxt(os.path.join(testdatadir, "slalib_hip_map.txt"))
rv = np.zeros_like(hip_tab['px'])

# Create V6C array from catalog data.
v6l = convert.cat2v6(hip_tab['raj2'], hip_tab['decj2'], hip_tab['pma'],
hip_tab['pmd'], hip_tab['px'], rv, tpm.CJ)

# UTC and TDB for mid-night of 2010/1/1.
utc = tpm.gcal2j(2010, 1, 1) - 0.5 # midnight
tt = tpm.utc2tdb(utc)

# Apply proper motion from J2000 to date.
v6o = convert.proper_motion(v6l, tt, tpm.J2000)

# Convert from mean equinox J2000 to true equinox and epoch of date.
v6o = convert.convertv6(v6o, s1=6, s2=11, utc=utc)

# Convert to Numpy rec-array.
cat = convert.v62cat(v6o, tpm.CJ)
cat = cat2array(cat)

ra = [math.degrees(j['alpha']) for j in cat]
dec = [math.degrees(j['delta']) for j in cat]
ra = np.array(ra)
dec = np.array(dec)

ra_diff = np.abs(ra - tab[:, 0]) * 3600.0
dec_diff = np.abs(dec - tab[:, 1]) * 3600.0
ra_diff = np.degrees(cat['alpha']) - tab[:, 0]
ra_diff = np.abs(ra_diff) * 3600.0
dec_diff = np.degrees(cat['delta']) - tab[:, 1]
dec_diff = np.abs(dec_diff) * 3600.0

print("Comparison with SLALIB map using HIPPARCOS data.")
fs = "{0} {1}\n" + \
"Min: {2:.4f} Max: {3:.4f} \nMean: {4:.4f} Std: {5:.4f}\n"
x = stats.describe(ra_diff)
Expand Down
15 changes: 8 additions & 7 deletions data/read_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -105,10 +105,11 @@ def __iter__(self):


def get_sla(filename):
f = open(os.path.join(testdatadir, filename), "r")
tab = csv.reader(CommentedFile(f), quoting=csv.QUOTE_NONNUMERIC,
delimiter=" ", skipinitialspace=True)

l = list(tab)
f.close()
return l
"""SLALIB ouput can have different columns."""
f = open(os.path.join(testdatadir, filename), "r")
tab = csv.reader(CommentedFile(f), quoting=csv.QUOTE_NONNUMERIC,
delimiter=" ", skipinitialspace=True)

l = np.array(list(tab), dtype=np.float64)
f.close()
return l
8 changes: 4 additions & 4 deletions data/slalib_convert_hip.py
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@
# for i, j in zip(ecl_lon, ecl_lat):
# # The format is very generous. The data are never this
# # accurate. .9 => ~1e-5 arc-sec.
# s = "%14.9f %14.9f \n"
# s = "%14.9f %14.9f\n"
# f.write(s % (i, j))

#sla_ecleq.
Expand All @@ -147,7 +147,7 @@
# for i, j in zip(raj2, decj2):
# # The format is very generous. The data are never this
# # accurate. .9 => ~1e-5 arc-sec.
# s = "%14.9f %14.9f \n"
# s = "%14.9f %14.9f\n"
# f.write(s % (i, j))

# sla_eqgal.
Expand All @@ -168,7 +168,7 @@
# for l, b in zip(gal_lon, gal_lat):
# # The format is very generous. The data are never this
# # accurate. .9 => ~1e-5 arc-sec.
# s = "%14.9f %14.9f \n"
# s = "%14.9f %14.9f\n"
# f.write(s % (l, b))

#sla_galeq.
Expand All @@ -189,7 +189,7 @@
# for r, d in zip(raj2, decj2):
# # The format is very generous. The data are never this
# # accurate. .9 => ~1e-5 arc-sec.
# s = "%14.9f %14.9f \n"
# s = "%14.9f %14.9f\n"
# f.write(s % (r, d))

# sla_map
Expand Down
Loading

0 comments on commit 4b4c4a7

Please sign in to comment.