Justine Devin / gammalib

Commit 1719d35f1a51248e4d7b68753f1a41dbcd4ff22a

Authored by Jürgen Knödlseder
2 parents 6f23e69f 1248ae77

Merge branch 'pep8' of https://cta-git.irap.omp.eu/gammalib into pep8 

Conflicts:
	test/test_GSky.py
Inline Side-by-side
Showing 1 changed file with 231 additions and 230 deletions
test/test_GSky.py
  View file @1719d35
... ... @@ -26,8 +26,8 @@ import os
26 26 # Text parser #
27 27 # =========== #
28 28 def loadtxt(filename):
29   - """
30   - Parse text file with columns of floats.
  29 + """
  30 + Parse text file with columns of floats.
31 31  
32 32 We have files with multiple spaces as delimiter,
33 33 which I think csv.reader can't handle.
... ... @@ -36,7 +36,7 @@ def loadtxt(filename):
36 36 data = []
37 37  
38 38 # Loop over all rows in file
39   - for row in [_.split() for _ in file(filename)]:
  39 + for row in [_.split() for _ in open(filename)]:
40 40 row_data = []
41 41 for item in row:
42 42 row_data.append(float(item))
... ... @@ -50,230 +50,231 @@ def loadtxt(filename):
50 50 # Test class for GammaLib sky module #
51 51 # ================================== #
52 52 class Test(GPythonTestSuite):
53   - """
54   - Test class for GammaLib sky module.
55   - """
56   - def __init__(self):
57   - """
58   - Constructor.
59   - """
60   - # Call base class constructor
61   - GPythonTestSuite.__init__(self)
62   -
63   - # Return
64   - return
65   -
66   - def set(self):
67   - """
68   - Set all test functions.
69   - """
70   - # Set test name
71   - self.name("sky")
72   -
73   - # Append tests
74   - self.append(self.test_skymap_healpix, "Test HEALPix map")
75   - self.append(self.test_skymap_car, "Test CAR projection map")
76   - self.append(self.test_skymap_tan, "Test TAN projection map")
77   - self.append(self.test_skymap_stg, "Test STG projection map")
78   - self.append(self.test_skymap_azp, "Test AZP projection map")
79   - self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion")
80   - # Return
81   - return
82   -
83   - # Generic skymap pixel transformation test
84   - def test_skymap_pixels(self, pixels, map):
85   - """
86   - Control coordinate and pixel transformation methods.
87   - """
88   - # Control pix2dir / dir2pix methods for all pixels
89   - for i in range(pixels.npix()):
90   - dir = pixels.pix2dir(i)
91   - ipix = pixels.dir2pix(dir)
92   - msg = map + " pix2dir/dir2pix check for pixel "+str(i)
93   - err = map + " GSkymap trouble with pixel "+str(i)+" ("+str(ipix)+ \
94   - "), RA="+str(dir.ra()*180/pi)+", Dec="+str(dir.dec()*180/pi)
95   - self.test_assert(i == ipix, msg, err)
96   -
97   - # Control SkyDir coordinate transformation for all pixels
98   - for i in range(pixels.npix()):
99   - dir = pixels.pix2dir(i)
100   - dir_new = GSkyDir()
101   - dir_new.lb(dir.l(),dir.b())
102   - dra = abs(dir.ra() - dir_new.ra())
103   - if (dra >= 5.0):
104   - dra -= 2.0*pi
105   - ddec = abs(dir.dec() - dir_new.dec())
106   - msg = map + " dir check for pixel "+str(i)
107   - err = map + " GSkymap trouble with pixel "+str(i)+" ("+str(dra)+ \
108   - ","+str(ddec)+")"
109   - self.test_assert(not (dra > 1.0e-9 or ddec > 1.0e-9), msg, err)
110   -
111   - # Return
112   - return
113   -
114   - # Generic skymap projection test
115   - def test_skymap_proj(self, proj):
116   - """
117   - Test skymap projection.
118   - """
119   - # Set filename
120   - file = "test_python_skymap_"+proj+".fits"
121   -
122   - # Remove test files
123   - try:
124   - os.remove(file)
125   - except:
126   - pass
127   -
128   - # Create skymap
129   - pixels = GSkymap(proj, "CEL", 83.6331, 22.0145, -3.7, 2.6, 5, 5, 20)
130   - for map in range(pixels.nmaps()):
131   - for i in range(pixels.npix()):
132   - pixels[i,map] = i+map*pixels.npix()
133   - pixels.save(file)
134   -
135   - # Test coordinate and pixel transformations
136   - self.test_skymap_pixels(pixels, proj)
137   -
138   - # Return
139   - return
140   -
141   - # Test HEALPix projection
142   - def test_skymap_healpix(self):
143   - """
144   - Test HEALPix interface.
145   - """
146   - # Set filenames
147   - file1 = "test_python_skymap_hpx_v1.fits"
148   - file2 = "test_python_skymap_hpx_v2.fits"
149   -
150   - # Remove test files
151   - try:
152   - os.remove(file1)
153   - os.remove(file2)
154   - except:
155   - pass
156   -
157   - # Create HEALPix skymap
158   - pixels = GSkymap("HPX", "GAL", 2, "RING", 2)
159   - for i in range(pixels.npix()):
160   - pixels[i] = i+1.0
161   - pixels[i,1] = i+1.0 + 1000.0
162   - pixels.save(file1)
163   -
164   - # Load HEALPix skymap
165   - pixels = GSkymap(file1)
166   -
167   - # Control coordinate and pixel transformations
168   - self.test_skymap_pixels(pixels, "HEALPix")
169   -
170   - # Save HEALPix skymap twice. The second saving should fail.
171   - try:
172   - pixels.save(file2, True)
173   - pixels.save(file2)
174   - except RuntimeError:
175   - pass
176   - else:
177   - raise RuntimeError("*** TEST ERROR: FITS file overwritten!")
178   - pixels.save(file2, True)
179   -
180   - # Load again HEALPix skymap
181   - pixels = GSkymap(file1)
182   -
183   - # Return
184   - return
185   -
186   - # Test CAR projection
187   - def test_skymap_car(self):
188   - """
189   - Test CAR projection.
190   - """
191   - # Execute generic test
192   - self.test_skymap_proj("CAR")
193   -
194   - # Return
195   - return
196   -
197   - # Test TAN projection
198   - def test_skymap_tan(self):
199   - """
200   - Test TAN projection.
201   - """
202   - # Execute generic test
203   - self.test_skymap_proj("TAN")
204   -
205   - # Return
206   - return
207   -
208   - # Test STG projection
209   - def test_skymap_stg(self):
210   - """
211   - Test STG projection.
212   - """
213   - # Execute generic test
214   - self.test_skymap_proj("STG")
215   -
216   - # Return
217   - return
218   -
219   - # Test AZP projection
220   - def test_skymap_azp(self):
221   - """
222   - Test AZP projection.
223   - """
224   - # Execute generic test
225   - self.test_skymap_proj("AZP")
226   -
227   - # Return
228   - return
229   -
230   - # FK5 J2000 to Galactic coordinate tranformations
231   - def test_fk5_to_galactic(self):
232   - """
233   - Test precision of FK5 J2000 to Galactic coordinate tranformations
234   - against pyast (http://dsberry.github.com/starlink/pyast.html)
235   - for 1000 random positions on the sky.
236   -
237   - Test txt files taken from:
238   - See https://github.com/astropy/coordinates-benchmark
239   - """
240   - # Set parameters
241   - fk5_filename = "data/initial_coords.txt"
242   - galactic_filename = "data/fk5_to_galactic.txt"
243   - fk5_coordinates = loadtxt(fk5_filename)
244   - galactic_coordinates = loadtxt(galactic_filename)
245   - coordinates = zip(fk5_coordinates, galactic_coordinates)
246   - max_dist = 0
247   -
248   - # Loop over coordinates
249   - for (ra, dec), (glon, glat) in coordinates:
250   -
251   - # Create point A in FK5 J2000 coordinates
252   - input_point = GSkyDir()
253   - input_point.radec_deg(ra, dec)
254   -
255   - # Convert point A to Galactic coordinates
256   - ll, bb = input_point.l_deg(), input_point.b_deg()
257   - actual_point = GSkyDir()
258   - actual_point.lb_deg(ll, bb)
259   -
260   - # Compute offset to reference result
261   - reference_point = GSkyDir()
262   - reference_point.lb_deg(glon, glat)
263   -
264   - # Compute distance and remember maximum
265   - this_dist = actual_point.dist_deg(reference_point)
266   - max_dist = max(max_dist, this_dist)
267   -
268   - # Convert max_dist from deg to milli-arcsec
269   - max_dist *= 1e3 * 3600
270   -
271   - # Results match within 10 milli-arcsec on my machine,
272   - # I put 20 arcsec here to avoid test failures on machines with
273   - # slightly different numerical results
274   - msg = ''
275   - err = ''
276   - self.test_assert(max_dist < 20, msg, err)
277   -
278   - # Return
279   - return
  53 + """
  54 + Test class for GammaLib sky module.
  55 + """
  56 + def __init__(self):
  57 + """
  58 + Constructor.
  59 + """
  60 + # Call base class constructor
  61 + GPythonTestSuite.__init__(self)
  62 +
  63 + # Return
  64 + return
  65 +
  66 + def set(self):
  67 + """
  68 + Set all test functions.
  69 + """
  70 + # Set test name
  71 + self.name("sky")
  72 +
  73 + # Append tests
  74 + self.append(self.test_skymap_healpix, "Test HEALPix map")
  75 + self.append(self.test_skymap_car, "Test CAR projection map")
  76 + self.append(self.test_skymap_tan, "Test TAN projection map")
  77 + self.append(self.test_skymap_stg, "Test STG projection map")
  78 + self.append(self.test_skymap_azp, "Test AZP projection map")
  79 + self.append(self.test_fk5_to_galactic, "Test FK5 to Galactic coordinate conversion")
  80 +
  81 + # Return
  82 + return
  83 +
  84 + # Generic skymap pixel transformation test
  85 + def test_skymap_pixels(self, pixels, map):
  86 + """
  87 + Control coordinate and pixel transformation methods.
  88 + """
  89 + # Control pix2dir / dir2pix methods for all pixels
  90 + for i in range(pixels.npix()):
  91 + dir = pixels.pix2dir(i)
  92 + ipix = pixels.dir2pix(dir)
  93 + msg = map + " pix2dir/dir2pix check for pixel " + str(i)
  94 + err = map + " GSkymap trouble with pixel " + str(i) + " (" + str(ipix) + \
  95 + "), RA=" + str(dir.ra() * 180 / pi) + ", Dec=" + str(dir.dec() * 180 / pi)
  96 + self.test_assert(i == ipix, msg, err)
  97 +
  98 + # Control SkyDir coordinate transformation for all pixels
  99 + for i in range(pixels.npix()):
  100 + dir = pixels.pix2dir(i)
  101 + dir_new = GSkyDir()
  102 + dir_new.lb(dir.l(), dir.b())
  103 + dra = abs(dir.ra() - dir_new.ra())
  104 + if (dra >= 5.0):
  105 + dra -= 2.0 * pi
  106 + ddec = abs(dir.dec() - dir_new.dec())
  107 + msg = map + " dir check for pixel " + str(i)
  108 + err = map + " GSkymap trouble with pixel " + str(i) + " (" + str(dra) + \
  109 + "," + str(ddec) + ")"
  110 + self.test_assert(not (dra > 1.0e-9 or ddec > 1.0e-9), msg, err)
  111 +
  112 + # Return
  113 + return
  114 +
  115 + # Generic skymap projection test
  116 + def test_skymap_proj(self, proj):
  117 + """
  118 + Test skymap projection.
  119 + """
  120 + # Set filename
  121 + file = "test_python_skymap_" + proj + ".fits"
  122 +
  123 + # Remove test files
  124 + try:
  125 + os.remove(file)
  126 + except:
  127 + pass
  128 +
  129 + # Create skymap
  130 + pixels = GSkymap(proj, "CEL", 83.6331, 22.0145, -3.7, 2.6, 5, 5, 20)
  131 + for map in range(pixels.nmaps()):
  132 + for i in range(pixels.npix()):
  133 + pixels[i, map] = i + map * pixels.npix()
  134 + pixels.save(file)
  135 +
  136 + # Test coordinate and pixel transformations
  137 + self.test_skymap_pixels(pixels, proj)
  138 +
  139 + # Return
  140 + return
  141 +
  142 + # Test HEALPix projection
  143 + def test_skymap_healpix(self):
  144 + """
  145 + Test HEALPix interface.
  146 + """
  147 + # Set filenames
  148 + file1 = "test_python_skymap_hpx_v1.fits"
  149 + file2 = "test_python_skymap_hpx_v2.fits"
  150 +
  151 + # Remove test files
  152 + try:
  153 + os.remove(file1)
  154 + os.remove(file2)
  155 + except:
  156 + pass
  157 +
  158 + # Create HEALPix skymap
  159 + pixels = GSkymap("HPX", "GAL", 2, "RING", 2)
  160 + for i in range(pixels.npix()):
  161 + pixels[i] = i + 1.0
  162 + pixels[i, 1] = i + 1.0 + 1000.0
  163 + pixels.save(file1)
  164 +
  165 + # Load HEALPix skymap
  166 + pixels = GSkymap(file1)
  167 +
  168 + # Control coordinate and pixel transformations
  169 + self.test_skymap_pixels(pixels, "HEALPix")
  170 +
  171 + # Save HEALPix skymap twice. The second saving should fail.
  172 + try:
  173 + pixels.save(file2, True)
  174 + pixels.save(file2)
  175 + except RuntimeError:
  176 + pass
  177 + else:
  178 + raise RuntimeError("*** TEST ERROR: FITS file overwritten!")
  179 + pixels.save(file2, True)
  180 +
  181 + # Load again HEALPix skymap
  182 + pixels = GSkymap(file1)
  183 +
  184 + # Return
  185 + return
  186 +
  187 + # Test CAR projection
  188 + def test_skymap_car(self):
  189 + """
  190 + Test CAR projection.
  191 + """
  192 + # Execute generic test
  193 + self.test_skymap_proj("CAR")
  194 +
  195 + # Return
  196 + return
  197 +
  198 + # Test TAN projection
  199 + def test_skymap_tan(self):
  200 + """
  201 + Test TAN projection.
  202 + """
  203 + # Execute generic test
  204 + self.test_skymap_proj("TAN")
  205 +
  206 + # Return
  207 + return
  208 +
  209 + # Test STG projection
  210 + def test_skymap_stg(self):
  211 + """
  212 + Test STG projection.
  213 + """
  214 + # Execute generic test
  215 + self.test_skymap_proj("STG")
  216 +
  217 + # Return
  218 + return
  219 +
  220 + # Test AZP projection
  221 + def test_skymap_azp(self):
  222 + """
  223 + Test AZP projection.
  224 + """
  225 + # Execute generic test
  226 + self.test_skymap_proj("AZP")
  227 +
  228 + # Return
  229 + return
  230 +
  231 + # FK5 J2000 to Galactic coordinate tranformations
  232 + def test_fk5_to_galactic(self):
  233 + """
  234 + Test precision of FK5 J2000 to Galactic coordinate tranformations
  235 + against pyast (http://dsberry.github.com/starlink/pyast.html)
  236 + for 1000 random positions on the sky.
  237 +
  238 + Test txt files taken from:
  239 + See https://github.com/astropy/coordinates-benchmark
  240 + """
  241 + # Set parameters
  242 + fk5_filename = "data/initial_coords.txt"
  243 + galactic_filename = "data/fk5_to_galactic.txt"
  244 + fk5_coordinates = loadtxt(fk5_filename)
  245 + galactic_coordinates = loadtxt(galactic_filename)
  246 + coordinates = zip(fk5_coordinates, galactic_coordinates)
  247 + max_dist = 0
  248 +
  249 + # Loop over coordinates
  250 + for (ra, dec), (glon, glat) in coordinates:
  251 +
  252 + # Create point A in FK5 J2000 coordinates
  253 + input_point = GSkyDir()
  254 + input_point.radec_deg(ra, dec)
  255 +
  256 + # Convert point A to Galactic coordinates
  257 + ll, bb = input_point.l_deg(), input_point.b_deg()
  258 + actual_point = GSkyDir()
  259 + actual_point.lb_deg(ll, bb)
  260 +
  261 + # Compute offset to reference result
  262 + reference_point = GSkyDir()
  263 + reference_point.lb_deg(glon, glat)
  264 +
  265 + # Compute distance and remember maximum
  266 + this_dist = actual_point.dist_deg(reference_point)
  267 + max_dist = max(max_dist, this_dist)
  268 +
  269 + # Convert max_dist from deg to milli-arcsec
  270 + max_dist *= 1e3 * 3600
  271 +
  272 + # Results match within 10 milli-arcsec on my machine,
  273 + # I put 20 arcsec here to avoid test failures on machines with
  274 + # slightly different numerical results
  275 + msg = ''
  276 + err = ''
  277 + self.test_assert(max_dist < 20, msg, err)
  278 +
  279 + # Return
  280 + return
... ...