我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用astropy.units.deg()。
def do_annual_parallax_test(filename): """testing functions called by a few unit tests""" with open(filename) as data_file: lines = data_file.readlines() ulens_params = lines[3].split() event_params = lines[4].split() data = np.loadtxt(filename, dtype=None) model = Model({ 't_0':float(ulens_params[1])+2450000., 'u_0':float(ulens_params[3]), 't_E':float(ulens_params[4]), 'pi_E_N':float(ulens_params[5]), 'pi_E_E':float(ulens_params[6]) }, coords=SkyCoord( event_params[1]+' '+event_params[2], unit=(u.deg, u.deg))) model.parameters.t_0_par = float(ulens_params[2])+2450000. time = data[:,0] dataset = MulensData([time, 20.+time*0., 0.1+time*0.,], add_2450000=True) model.set_datasets([dataset]) model.parallax(satellite=False, earth_orbital=True, topocentric=False) return np.testing.assert_almost_equal( model.data_magnification[0] / data[:,1], 1.0, decimal=4)
def main(): parser = argparse.ArgumentParser( description="Convert among multiple coordinate formats") parser.add_argument("coord", nargs="+") args = parser.parse_args() ra, dec = parse_coord(args.coord) info = ( "%-14s %-14s\n" % ("R.A.", "Dec.") + "%s--%s\n" % ("-"*14, "-"*14) + "%-14.3f %-+14.3f\n" % (ra.deg, dec.deg) + "%-14s %-14s\n" % ( ra.to_string(unit=au.hourangle, precision=4), dec.to_string(unit=au.deg, alwayssign=True, precision=3)) + "%-14s %-14s\n" % ( ra.to_string(unit=au.hourangle, sep=":", precision=4), dec.to_string(unit=au.deg, alwayssign=True, sep=":", precision=3)) + "%-14s %-14s\n" % ( ra.to_string(unit=au.hourangle, sep=" ", precision=4), dec.to_string(unit=au.deg, alwayssign=True, sep=" ", precision=3)) ) print(info)
def test_meta_from_position(igmsp): # One match meta = igmsp.meta_from_position((0.0019,17.7737), 1*u.arcsec) assert len(meta) == 1 # Blank meta2 = igmsp.meta_from_position((10.038604,55.298477), 1*u.arcsec) assert meta2 is None # Multiple sources (insure rank order) meta3 = igmsp.meta_from_position((0.0055,-1.5), 1*u.deg) assert len(meta3) == 2 assert np.isclose(meta3['R'][0],meta3['R'][1]) # Multiple meta entries (GGG) meta4 = igmsp.meta_from_position('001115.23+144601.8', 1*u.arcsec) assert len(meta4) == 2 assert meta4['R'][0] != meta4['R'][1] # Multiple but grab closest source meta5 = igmsp.meta_from_position((0.0055,-1.5), 1*u.deg, max_match=1) assert len(meta5) == 1 # Groups meta = igmsp.meta_from_position((2.813500,14.767200), 20*u.deg, groups=['GGG','HD-LLS_DR1']) for group in meta['GROUP'].data: assert group in ['GGG', 'HD-LLS_DR1'] # Physical separation meta6 = igmsp.meta_from_position('001115.23+144601.8', 300*u.kpc) assert len(meta6) == 2
def test_query_position(igmsp): # One match _, _, idx = igmsp.qcat.query_position((0.0019,17.7737), 1*u.arcsec) assert idx >= 0 # Blank _, _, idx = igmsp.qcat.query_position((10.038604,55.298477), 1*u.arcsec) assert len(idx) == 0 # Multiple (insure rank order) icoord = SkyCoord(ra=0.0055, dec=-1.5, unit='deg') _, subcat, _ = igmsp.qcat.query_position(icoord, 1*u.deg) # Test coord = SkyCoord(ra=subcat['RA'], dec=subcat['DEC'], unit='deg') sep = icoord.separation(coord) isrt = np.argsort(sep) assert isrt[0] == 0 assert isrt[-1] == len(subcat)-1 # Multiple but grab only 1 _, _, idxs = igmsp.qcat.query_position(icoord, 1*u.deg, max_match=1) assert len(idxs) == 1
def test_cat_from_query_coord(igmsp): # Single coord = SkyCoord(ra=0.0019, dec=17.7737, unit='deg') # _, ccat, _ = igmsp.qcat.query_coords(coord) assert ccat['IGM_ID'][0] == 0 # Multiple coords = SkyCoord(ra=[0.0028,0.0019], dec=[14.9747,17.7737], unit='deg') _, ccat2, _ = igmsp.qcat.query_coords(coords) assert len(ccat2) == 2 assert ccat2['IGM_ID'][0] == 1 assert ccat2['IGM_ID'][1] == 0 # One miss coords3 = SkyCoord(ra=[9.0028,0.0019], dec=[-14.9747,17.7737], unit='deg') _, ccat3, _ = igmsp.qcat.query_coords(coords3) assert ccat3['IGM_ID'][0] == -1
def test_spectra_from_coord(igmsp): # No match specN, metaN = igmsp.spectra_from_coord((0.0019, -17.7737)) assert specN is None assert metaN is None # One match spec, meta = igmsp.spectra_from_coord((0.0019, 17.7737)) assert spec.nspec == 1 assert meta['PLATE'][0] == 6173 # Multiple matches and spectra spec, meta = igmsp.spectra_from_coord('001115.23+144601.8')#, groups=['GGG']) # Not in debug file for BOSS or SDSS assert spec.nspec == 2 # Many matches; takes closest spec, meta = igmsp.spectra_from_coord((0.0019, 17.7737), tol=10*u.deg) assert spec.nspec == 1 assert meta['PLATE'][0] == 6173
def setUp(self): self.dir = './test_results' self.lowcore = create_named_configuration('LOWBD2-CORE') os.makedirs(self.dir, exist_ok=True) self.times = (numpy.pi / (12.0)) * numpy.linspace(-3.0, 3.0, 7) self.frequency = numpy.array([1e8]) self.channel_bandwidth = numpy.array([1e6]) self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000') self.vis = create_visibility(self.lowcore, self.times, self.frequency, channel_bandwidth=self.channel_bandwidth, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame('stokesI')) self.vis.data['vis'] *= 0.0 # Create model self.model = create_test_image(cellsize=0.0015, phasecentre=self.vis.phasecentre, frequency=self.frequency) self.model.data[self.model.data > 1.0] = 1.0 self.vis = predict_2d(self.vis, self.model) assert numpy.max(numpy.abs(self.vis.vis)) > 0.0 export_image_to_fits(self.model, '%s/test_solve_skycomponent_model.fits' % (self.dir)) self.bigmodel = create_image_from_visibility(self.vis, cellsize=0.0015, npixel=512)
def _checkcomponents(self, dirty, fluxthreshold=5.0, positionthreshold=1.0): comps = find_skycomponents(dirty, fwhm=1.0, threshold=fluxthreshold, npixels=5) assert len(comps) == len(self.components), "Different number of components found: original %d, recovered %d" % \ (len(self.components), len(comps)) cellsize = abs(dirty.wcs.wcs.cdelt[0]) # Check for agreement between image and DFT for comp in comps: sflux = sum_visibility(self.componentvis, comp.direction)[0] assert abs(comp.flux[0, 0] - sflux[0, 0]) < fluxthreshold, \ "Fitted and DFT flux differ %s %s" % (comp.flux[0, 0], sflux[0, 0]) # Check for agreement in direction ocomp = find_nearest_component(comp.direction, self.components) radiff = abs(comp.direction.ra.deg - ocomp.direction.ra.deg) / cellsize assert radiff < positionthreshold, "Component differs in dec %.3f pixels" % radiff decdiff = abs(comp.direction.dec.deg - ocomp.direction.dec.deg) / cellsize assert decdiff < positionthreshold, "Component differs in dec %.3f pixels" % decdiff
def actualSetup(self, sky_pol_frame='stokesIQUV', data_pol_frame='linear', f=None): if f is None: f = [100.0, 50.0, -10.0, 40.0] f = numpy.array(f) self.flux = numpy.array([f, 0.8 * f, 0.6 * f]) # The phase centre is absolute and the component is specified relative (for now). # This means that the component should end up at the position phasecentre+compredirection self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000') self.compabsdirection = SkyCoord(ra=+181.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000') self.comp = Skycomponent(direction=self.compabsdirection, frequency=self.frequency, flux=self.flux, polarisation_frame=PolarisationFrame(sky_pol_frame)) self.vis = create_blockvisibility(self.lowcore, self.times, self.frequency, phasecentre=self.phasecentre, channel_bandwidth=self.channel_bandwidth, weight=1.0, polarisation_frame=PolarisationFrame(data_pol_frame)) self.vis = predict_skycomponent_blockvisibility(self.vis, self.comp)
def setUp(self): self.dir = './test_results' os.makedirs(self.dir, exist_ok=True) self.vnchan = 5 self.lowcore = create_named_configuration('LOWBD2-CORE') self.times = (numpy.pi / 12.0) * numpy.linspace(-3.0, 3.0, 7) self.frequency = numpy.linspace(8e7, 1.2e8, self.vnchan) self.startfrequency = numpy.array([8e7]) self.channel_bandwidth = numpy.array(self.vnchan * [self.frequency[1] - self.frequency[0]]) self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000') self.vis = create_visibility(self.lowcore, times=self.times, frequency=self.frequency, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame('stokesI'), channel_bandwidth=self.channel_bandwidth) self.model = create_image_from_visibility(self.vis, npixel=512, cellsize=0.001, nchan=self.vnchan, frequency=self.startfrequency)
def setUp(self): self.lowcore = create_named_configuration('LOWBD2-CORE') self.times = (numpy.pi / 43200.0) * numpy.arange(0.0, 300.0, 30.0) self.frequency = numpy.linspace(1.0e8, 1.1e8, 3) self.channel_bandwidth = numpy.array([1e7, 1e7, 1e7]) # Define the component and give it some spectral behaviour f = numpy.array([100.0, 20.0, -10.0, 1.0]) self.flux = numpy.array([f, 0.8 * f, 0.6 * f]) # The phase centre is absolute and the component is specified relative (for now). # This means that the component should end up at the position phasecentre+compredirection self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000') self.compabsdirection = SkyCoord(ra=+181.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000') pcof = self.phasecentre.skyoffset_frame() self.compreldirection = self.compabsdirection.transform_to(pcof) self.comp = Skycomponent(direction=self.compreldirection, frequency=self.frequency, flux=self.flux)
def test_phase_rotation_identity(self): self.vis = create_visibility(self.lowcore, self.times, self.frequency, channel_bandwidth=self.channel_bandwidth, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame("stokesIQUV")) self.vismodel = predict_skycomponent_visibility(self.vis, self.comp) newphasecenters = [SkyCoord(182, -35, unit=u.deg), SkyCoord(182, -30, unit=u.deg), SkyCoord(177, -30, unit=u.deg), SkyCoord(176, -35, unit=u.deg), SkyCoord(216, -35, unit=u.deg), SkyCoord(180, -70, unit=u.deg)] for newphasecentre in newphasecenters: # Phase rotating back should not make a difference original_vis = self.vismodel.vis original_uvw = self.vismodel.uvw rotatedvis = phaserotate_visibility(phaserotate_visibility(self.vismodel, newphasecentre, tangent=False), self.phasecentre, tangent=False) assert_allclose(rotatedvis.uvw, original_uvw, rtol=1e-7) assert_allclose(rotatedvis.vis, original_vis, rtol=1e-7)
def setUp(self): self.dir = './test_results' self.lowcore = create_named_configuration('LOWBD2-CORE') os.makedirs(self.dir, exist_ok=True) self.times = (numpy.pi / (12.0)) * numpy.linspace(-3.0, 3.0, 7) self.nchan=8 self.frequency = numpy.linspace(0.8e8, 1.2e8, self.nchan) self.channel_bandwidth = numpy.array(self.nchan*[self.frequency[1]-self.frequency[0]]) self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000') self.vis = create_visibility(self.lowcore, self.times, self.frequency, channel_bandwidth=self.channel_bandwidth, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame('stokesI')) self.vis.data['vis'] *= 0.0 # Create model self.model = create_test_image(cellsize=0.0015, phasecentre=self.vis.phasecentre, frequency=self.frequency) self.model.data[self.model.data > 1.0] = 1.0 self.vis = predict_2d(self.vis, self.model) assert numpy.max(numpy.abs(self.vis.vis)) > 0.0 export_image_to_fits(self.model, '%s/test_solve_image_mm_model.fits' % (self.dir)) self.bigmodel = create_image_from_visibility(self.vis, cellsize=0.0015, npixel=512, frequency=self.frequency)
def random_points(self, n=1): """ Generate uniformly distributed random points within the sky image (coverage). Parameters ---------- n : int, optional The number of random points required. Default: 1 Returns ------- lon, lat : float, or 1D `~numpy.ndarray` The longitudes and latitudes (in world coordinate) generated. Unit: [deg] """ raise NotImplementedError
def world2pix(self, x, y): """ Convert the world coordinates (R.A., Dec.) into the pixel coordinates (indexes) within the sky data array. Parameters ---------- x, y : float, `~numpy.ndarray` The R.A., Dec. world coordinates Unit: [deg] Returns ------- ri, ci : int, `~numpy.ndarray` The row, column indexes within the sky data array. """ pixelsize = self.pixelsize * AUC.arcsec2deg # [deg] x, y = np.asarray(x), np.asarray(y) # [deg] ri0, ci0 = self.ysize//2, self.xsize//2 ri = np.round((y - self.ycenter) / pixelsize + ri0).astype(int) ci = np.round((x - self.xcenter) / pixelsize + ci0).astype(int) return (ri, ci)
def random_points(self, n=1): """ Generate uniformly distributed random points within the sky patch. Returns ------- lon : float, or 1D `~numpy.ndarray` Longitudes (Galactic/equatorial); Unit: [deg] lat : float, or 1D `~numpy.ndarray` Latitudes (Galactic/equatorial); Unit: [deg] """ lon_min, lon_max = self.lon_limit lat_min, lat_max = self.lat_limit lon = np.random.uniform(low=lon_min, high=lon_max, size=n) lat = np.random.uniform(low=lat_min, high=lat_max, size=n) return (lon, lat)
def signal_gaussian( signal_location=np.array([61, 21])*u.deg, fov_center=np.array([60, 20])*u.deg, width=0.05*u.deg, signal_events=80, bins=[80, 80], fov=4*u.deg, ): # reshape so if signal_events = 1 the array can be indexed in the same way. signal = multivariate_normal.rvs( mean=signal_location.value, cov=width.value, size=signal_events ).reshape(signal_events, 2) r = np.array([fov_center - fov/2, fov_center + fov/2]).T signal_hist, _, _ = np.histogram2d(signal[:, 0], signal[:, 1], bins=bins, range=r) return signal_hist
def prepare(self, path, star_finder): with warnings.catch_warnings(): warnings.simplefilter('ignore', AstropyWarning) self.calibration.selectImage(path) self.names, self.vmag, alt, az = self.calibration.catalog.filter(Configuration.min_alt * u.deg, Configuration.max_mag) altaz = np.array([alt.radian, az.radian]).transpose() self.pos = np.column_stack(self.calibration.project(altaz)) self.finder = star_finder self.image = cv2.imread(path) self.finder.setImage(self.image) self.altaz = np.array([alt.degree, az.degree]).transpose()
def galToCel(ll, bb): """ Converts Galactic (deg) to Celestial J2000 (deg) coordinates """ bb = numpy.radians(bb) sin_bb = numpy.sin(bb) cos_bb = numpy.cos(bb) ll = numpy.radians(ll) ra_gp = numpy.radians(192.85948) de_gp = numpy.radians(27.12825) lcp = numpy.radians(122.932) sin_lcp_ll = numpy.sin(lcp - ll) cos_lcp_ll = numpy.cos(lcp - ll) sin_d = (numpy.sin(de_gp) * sin_bb) \ + (numpy.cos(de_gp) * cos_bb * cos_lcp_ll) ramragp = numpy.arctan2(cos_bb * sin_lcp_ll, (numpy.cos(de_gp) * sin_bb) \ - (numpy.sin(de_gp) * cos_bb * cos_lcp_ll)) dec = numpy.arcsin(sin_d) ra = (ramragp + ra_gp + (2. * numpy.pi)) % (2. * numpy.pi) return numpy.degrees(ra), numpy.degrees(dec)
def dec2dms(dec): """ ADW: This should really be replaced by astropy """ DEGREE = 360. HOUR = 24. MINUTE = 60. SECOND = 3600. if isinstance(dec,basestring): dec = float(dec) sign = numpy.copysign(1.0,dec) fdeg = np.abs(dec) deg = int(fdeg) fminute = (fdeg - deg)*MINUTE minute = int(fminute) second = (fminute - minute)*MINUTE deg = int(deg * sign) return (deg, minute, second)
def test_xarray(): import numpy as np import astropy.units as u from ..utils import xarray, yarray x = xarray((3, 3)) xx = np.array([[0, 1, 2], [0, 1, 2], [0, 1, 2]]) assert np.all(x == xx) x = xarray((3, 3), yx=(1, 0)) assert np.all(x == xx) x = xarray((3, 3), yx=(0, 0.5)) xx = xx - 0.5 assert np.all(x == xx) x = xarray((3, 3), th=90 * u.deg) y = yarray((3, 3)) assert np.allclose(x, y)
def test_yarray(): import numpy as np import astropy.units as u from ..utils import yarray, xarray y = yarray((3, 3)) yy = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]]) assert np.all(y == yy) y = yarray((3, 3), yx=(0, 1)) assert np.all(y == yy) y = yarray((3, 3), yx=(0.5, 0)) yy = yy - 0.5 assert np.all(y == yy) y = yarray((3, 3), th=-90 * u.deg) x = xarray((3, 3)) assert np.allclose(y, x)
def match_to_database(source_list, image_list, cat, Ast, imscale, thresh=8.): ''' RA, Decdetection catalogs and vizier (astroquery) match? vo conesearch? ''' matches = [] for sources, image in zip(source_list, image_list): ifuslot = sources[1] Ast.get_ifuslot_projection(ifuslot, imscale, image.shape[1]/2., image.shape[0]/2.) for source in sources[0]: ra, dec = Ast.tp_ifuslot.wcs.wcs_pix2world(source['xcentroid'], source['ycentroid'],1) c = SkyCoord(ra, dec, unit=(unit.degree, unit.degree)) idx, d2d, d3d = match_coordinates_sky(c, cat, nthneighbor=1) if d2d.arcsec[0] < thresh: dra = cat.ra[idx].deg - float(ra) ddec = cat.dec[idx].deg - float(dec) matches.append([int(ifuslot),int(idx), float(ra), float(dec), cat.ra[idx].deg, cat.dec[idx].deg, dra, ddec, d2d.arcsec[0]]) return matches
def get_wise(ra,dec,band): mission='wise' dataset='allwise' table='p3am_cdd' successful=False while not successful: try: results=Ibe.query_region(coord.SkyCoord(ra, dec, unit=(u.deg, u.deg), frame='icrs'),mission=mission,dataset=dataset,table=table) successful=True except requests.exceptions.ConnectionError: print 'Connection failed, retrying' sleep(10) url = 'http://irsa.ipac.caltech.edu/ibe/data/'+mission+'/'+dataset+'/'+table+'/' params = { 'coadd_id': results[results['band']==band]['coadd_id'][0], 'band': band } params['coaddgrp'] = params['coadd_id'][:2] params['coadd_ra'] = params['coadd_id'][:4] path = str.format('{coaddgrp:s}/{coadd_ra:s}/{coadd_id:s}/{coadd_id:s}-w{band:1d}-int-3.fits',**params) outname=path.split('/')[-1] download_file(url+path,outname) return outname
def peek(self, figsize=(10, 10), color='b', alpha=0.75, print_to_file=None, **kwargs): """ Show extracted fieldlines overlaid on HMI image. """ fig = plt.figure(figsize=figsize) ax = fig.gca(projection=self.hmi_map) self.hmi_map.plot() ax.set_autoscale_on(False) for stream, _ in self.streamlines: ax.plot(self._convert_angle_to_length(stream[:, 0]*u.cm, working_units=u.arcsec).to(u.deg), self._convert_angle_to_length(stream[:, 1]*u.cm, working_units=u.arcsec).to(u.deg), alpha=alpha, color=color, transform=ax.get_transform('world')) if print_to_file is not None: plt.savefig(print_to_file, **kwargs) plt.show()
def coords(self): """Coordinates of event""" coords = SkyCoord( self.event_params[1]+' '+self.event_params[2], unit=(u.deg, u.deg)) return coords
def test_model_coords(): coords = SkyCoord('18:00:00 -30:00:00', unit=(u.hourangle, u.deg)) model_1 = Model({'t_0':2450000, 'u_0':0.1, 't_E':100}, coords='18:00:00 -30:00:00') assert isinstance(model_1.coords, SkyCoord) assert model_1.coords.ra == coords.ra assert model_1.coords.dec == coords.dec assert model_1.coords.dec.deg == -30.00 model_3 = Model({'t_0':2450000, 'u_0':0.1, 't_E':100}) model_3.coords = '17:00:00 -27:32:14' assert model_3.coords.to_string('hmsdms') == '17h00m00s -27d32m14s'
def test_data_coords(): coords = SkyCoord('18:00:00 -30:00:00', unit=(u.hourangle, u.deg)) data_1 = MulensData( file_name=SAMPLE_FILE_01, coords='18:00:00 -30:00:00') assert isinstance(data_1.coords, SkyCoord) assert data_1.coords.ra == coords.ra assert data_1.coords.dec == coords.dec assert data_1.coords.dec.deg == -30.00 data_3 = MulensData(file_name=SAMPLE_FILE_01) data_3.coords = '17:00:00 -27:32:14' assert data_3.coords.to_string('hmsdms') == '17h00m00s -27d32m14s'
def test_BLPS_01(): """simple binary lens with point source""" params = ModelParameters({ 't_0':6141.593, 'u_0':0.5425, 't_E':62.63*u.day, 'alpha':49.58*u.deg, 's':1.3500, 'q':0.00578}) model = Model(parameters=params) t = np.array([6112.5]) data = MulensData(data_list=[t, t*0.+16., t*0.+0.01]) model.set_datasets([data]) magnification = model.data_magnification[0][0] np.testing.assert_almost_equal(magnification, 4.691830781584699) # This value comes from early version of this code. # np.testing.assert_almost_equal(m, 4.710563917) # This value comes from Andy's getbinp().
def test_BLPS_02_AC(): """simple binary lens with extended source and different methods to evaluate magnification - version with adaptivecontouring """ params = ModelParameters({ 't_0':2456141.593, 'u_0':0.5425, 't_E':62.63*u.day, 'alpha':49.58*u.deg, 's':1.3500, 'q':0.00578, 'rho':0.01}) model = Model(parameters=params) t = (np.array([6112.5, 6113., 6114., 6115., 6116., 6117., 6118., 6119]) + 2450000.) ac_name = 'AdaptiveContouring' methods = [2456113.5, 'Quadrupole', 2456114.5, 'Hexadecapole', 2456116.5, ac_name, 2456117.5] accuracy_1 = {'accuracy': 0.04} accuracy_2 = {'accuracy': 0.01, 'ld_accuracy': 0.0001} model.set_magnification_methods(methods) model.set_magnification_methods_parameters({ac_name: accuracy_1}) data = MulensData(data_list=[t, t*0.+16., t*0.+0.01]) model.set_datasets([data]) result = model.data_magnification[0] expected = np.array([4.69183078, 2.87659723, 1.83733975, 1.63865704, 1.61038135, 1.63603122, 1.69045492, 1.77012807]) np.testing.assert_almost_equal(result, expected, decimal=3) # Below we test passing the limb coef to VBBL function. data.bandpass = 'I' model.set_limb_coeff_u('I', 10.) # This is an absurd value but I needed something quick. model.set_magnification_methods_parameters({ac_name: accuracy_2}) result = model.data_magnification[0] np.testing.assert_almost_equal(result[5], 1.6366862, decimal=3)
def test_methods_parameters(): """make sure additional parameters are properly passed to very inner functions""" params = ModelParameters({ 't_0':2456141.593, 'u_0':0.5425, 't_E':62.63*u.day, 'alpha':49.58*u.deg, 's':1.3500, 'q':0.00578, 'rho':0.01}) model = Model(parameters=params) t = np.array([2456117.]) methods = [2456113.5, 'Quadrupole', 2456114.5, 'Hexadecapole', 2456116.5, 'VBBL', 2456117.5] model.set_magnification_methods(methods) data = MulensData(data_list=[t, t*0.+16., t*0.+0.01]) model.set_datasets([data]) result_1 = model.data_magnification[0] vbbl_options = {'accuracy': 0.1} methods_parameters = {'VBBL': vbbl_options} model.set_magnification_methods_parameters(methods_parameters) result_2 = model.data_magnification[0] vbbl_options = {'accuracy': 1.e-5} methods_parameters = {'VBBL': vbbl_options} model.set_magnification_methods_parameters(methods_parameters) result_3 = model.data_magnification[0] assert result_1[0] != result_2[0] assert result_1[0] != result_3[0] assert result_2[0] != result_3[0]
def __init__(self, *args, **kwargs): if not isinstance(args[0], SkyCoord) and 'unit' not in kwargs: kwargs['unit'] = (u.hourangle, u.deg) SkyCoord.__init__(self, *args, **kwargs)
def galactic_l(self): """ Galactic longitude. Note that for connivance, the values l > 180 degrees are represented as 360-l. """ l = self.galactic.l if l > 180. * u.deg: l = l - 360. * u.deg return l
def dalpha_dt(self): """ change in alpha vs. time in deg/yr""" return self.parameters['dalpha_dt'].to(u.deg / u.yr).value
def xyz(self): """ *Astropy.CartesianRepresentation* return X,Y,Z positions based on RA, DEC and distance """ if self._xyz is None: ra_dec = [ text.decode('UTF-8') for text in self.data_lists['ra_dec']] self._xyz = SkyCoord( ra_dec, distance=self.data_lists['distance'], unit=(u.hourangle, u.deg, u.au)).cartesian return self._xyz
def parse_coord(c): if len(c) == 6: # h m s d m s ra = Angle((float(c[0]), float(c[1]), float(c[2])), unit=au.hourangle) dec = Angle((float(c[3]), float(c[4]), float(c[5])), unit=au.deg) elif len(c) == 2: ra = Angle(float(c[0]), unit=au.deg) dec = Angle(float(c[1]), unit=au.deg) else: raise ValueError("invalid coordinate: {0}".format(c)) return (ra, dec)
def input_params(votbl=None): """ Parameters ---------- votbl Returns ------- """ if votbl is None: votbl = empty_vo(rtype='meta') all_params = [] # POS pos = Param(votbl, name="INPUT:POS", value="", datatype="char", arraysize="*") pos.description = ('The center of the region of interest. The coordinate values are '+ 'specified in list format (comma separated) in decimal degrees with '+ 'no embedded white space followed by an optional coord. system. '+ 'Allowed systems are (ICRS) and the default is ICRS.') all_params.append(pos) # SIZE size = Param(votbl, name="INPUT:SIZE", value="0.1", datatype="double", unit="deg") size.description = ('The radius of the circular region of interest in decimal degrees.'+ 'Default sized radius is 0.001 degrees') all_params.append(size) # BAND band = Param(votbl, name="INPUT:BAND", value="ALL", datatype="char", arraysize="*") band.description = 'Not currently implemented' all_params.append(band) # TIME ptime = Param(votbl, name="INPUT:TIME", value="", datatype="char", arraysize="*") ptime.description = 'Not currently implemented' all_params.append(ptime) # FORMAT format = Param(votbl, name="INPUT:FORMAT", value="ALL", datatype="char", arraysize="*") format.description = ('Desired format of retrieved data. \n'+ 'Allowed values are HDF5, METADATA') all_params.append(format) # Return return all_params
def test_query_coords(igmsp): # Single coord = SkyCoord(ra=0.0019, dec=17.7737, unit='deg') # _, _, idx = igmsp.qcat.query_coords(coord) assert idx[0] >= 0 # Multiple coords = SkyCoord(ra=[0.0019]*2, dec=[17.7737]*2, unit='deg') _, _, idxs = igmsp.qcat.query_coords(coords) assert len(idxs) == 2 # Dataset _, _, idxs2 = igmsp.qcat.query_coords(coords, groups=['BOSS_DR12']) assert np.sum(idxs2 >= 0) == 2 _, _, idxs3 = igmsp.qcat.query_coords(coords, groups=['HD-LLS_DR1']) assert np.sum(idxs3 >= 0) == 0
def test_spectra_in_group(igmsp): # Missed a source -- raises IOError coords = SkyCoord(ra=[0.0028, 0.0019], dec=[14.9747, -17.77374], unit='deg') with pytest.raises(IOError): spec, meta = igmsp.spectra_in_group(coords, 'BOSS_DR12') # Another with both missing coords = SkyCoord(ra=[2.8135,16.5802], dec=[-14.7672, -0.8065], unit='deg') with pytest.raises(IOError): spec, meta = igmsp.spectra_in_group(coords, 'GGG') # Each source has only spectrum in the group coords = SkyCoord(ra=[0.0028, 0.0019], dec=[14.9747, 17.77374], unit='deg') spec, meta = igmsp.spectra_in_group(coords, 'BOSS_DR12') # Test assert spec.nspec == 2 assert meta['PLATE'][0] == 6177 # Each source has multiple spectra in the group coords = SkyCoord(ra=[2.8135,16.5802], dec=[14.7672, 0.8065], unit='deg') spec, meta = igmsp.spectra_in_group(coords, 'GGG') assert meta['DISPERSER'][0] == 'B600' qdict = dict(DISPERSER='R400') spec, meta = igmsp.spectra_in_group(coords, 'GGG', query_dict=qdict) assert meta['DISPERSER'][0] == 'R400' # Another with bad grating qdict = dict(DISPERSER='X400') with pytest.raises(IOError): spec, meta = igmsp.spectra_in_group(coords, 'GGG', query_dict=qdict) ''' # Multiple spectra per group coords = SkyCoord(ra=[2.8135, 16.5802], dec=[14.7672, 0.8065], unit='deg') spec, meta = igmsp.coords_to_spectra(coords, 'GGG', all_spec=True) assert spec.nspec == 4 '''
def setUp(self): self.dir = './test_results' os.makedirs(self.dir, exist_ok=True) self.frequency = numpy.linspace(1e8, 1.5e8, 3) self.channel_bandwidth = numpy.array([2.5e7, 2.5e7, 2.5e7]) self.flux = numpy.array([[100.0], [100.0], [100.0]]) self.phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000') self.config = create_named_configuration('LOWBD2-CORE') self.times = numpy.linspace(-300.0, 300.0, 3) * numpy.pi / 43200.0 nants = self.config.xyz.shape[0] assert nants > 1 assert len(self.config.names) == nants assert len(self.config.mount) == nants
def createVis(self, config, dec=-35.0, rmax=None): self.config = create_named_configuration(config, rmax=rmax) self.phasecentre = SkyCoord(ra=+15 * u.deg, dec=dec * u.deg, frame='icrs', equinox='J2000') self.vis = create_visibility(self.config, self.times, self.frequency, channel_bandwidth=self.channel_bandwidth, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame('stokesI'))
def test_insert_skycomponent_nearest(self): dphasecentre = SkyCoord(ra=+181.0 * u.deg, dec=-58.0 * u.deg, frame='icrs', equinox='J2000') sc = create_skycomponent(direction=dphasecentre, flux=numpy.array([[1.0]]), frequency=self.frequency, polarisation_frame=PolarisationFrame('stokesI')) self.model.data *= 0.0 insert_skycomponent(self.model, sc, insert_method='Nearest') # These test a regression but are not known a priori to be correct self.assertAlmostEqual(self.model.data[0, 0, 151, 122], 1.0, 7) self.assertAlmostEqual(self.model.data[0, 0, 152, 122], 0.0, 7)
def test_insert_skycomponent_sinc(self): dphasecentre = SkyCoord(ra=+181.0 * u.deg, dec=-58.0 * u.deg, frame='icrs', equinox='J2000') sc = create_skycomponent(direction=dphasecentre, flux=numpy.array([[1.0]]), frequency=self.frequency, polarisation_frame=PolarisationFrame('stokesI')) self.model.data *= 0.0 insert_skycomponent(self.model, sc, insert_method='Sinc') # These test a regression but are not known a priori to be correct self.assertAlmostEqual(self.model.data[0, 0, 151, 122], 0.87684398703184396, 7) self.assertAlmostEqual(self.model.data[0, 0, 152, 122], 0.2469311811046056, 7)
def test_insert_skycomponent_lanczos(self): dphasecentre = SkyCoord(ra=+181.0 * u.deg, dec=-58.0 * u.deg, frame='icrs', equinox='J2000') sc = create_skycomponent(direction=dphasecentre, flux=numpy.array([[1.0]]), frequency=self.frequency, polarisation_frame=PolarisationFrame('stokesI')) self.model.data *= 0.0 insert_skycomponent(self.model, sc, insert_method='Lanczos') # These test a regression but are not known a priori to be correct self.assertAlmostEqual(self.model.data[0,0,151, 122], 0.87781267543090036, 7) self.assertAlmostEqual(self.model.data[0,0,152, 122], 0.23817562762032077, 7)
def test_insert_skycomponent_lanczos_bandwidth(self): dphasecentre = SkyCoord(ra=+181.0 * u.deg, dec=-58.0 * u.deg, frame='icrs', equinox='J2000') sc = create_skycomponent(direction=dphasecentre, flux=numpy.array([[1.0]]), frequency=self.frequency, polarisation_frame=PolarisationFrame('stokesI')) self.model.data *= 0.0 insert_skycomponent(self.model, sc, insert_method='Lanczos', bandwidth=0.5) # These test a regression but are not known a priori to be correct self.assertAlmostEqual(self.model.data[0,0,151, 122], 0.24031092091707615, 7) self.assertAlmostEqual(self.model.data[0,0,152, 122], 0.18648989466050975, 7)
def setUp(self): self.dir = './test_results' os.makedirs(self.dir, exist_ok=True) self.lowcore = create_named_configuration('LOWBD2-CORE') self.times = numpy.linspace(-3, +3, 13) * (numpy.pi / 12.0) self.frequency = numpy.array([1e8]) self.channel_bandwidth = numpy.array([1e7]) # Define the component and give it some polarisation and spectral behaviour f = numpy.array([100.0]) self.flux = numpy.array([f]) self.phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000') self.compabsdirection = SkyCoord(ra=17.0 * u.deg, dec=-36.5 * u.deg, frame='icrs', equinox='J2000') self.comp = create_skycomponent(flux=self.flux, frequency=self.frequency, direction=self.compabsdirection, polarisation_frame=PolarisationFrame('stokesI')) self.image = create_test_image(frequency=self.frequency, phasecentre=self.phasecentre, cellsize=0.001, polarisation_frame=PolarisationFrame('stokesI')) self.image.data[self.image.data < 0.0] = 0.0 self.image_graph = delayed(create_test_image)(frequency=self.frequency, phasecentre=self.phasecentre, cellsize=0.001, polarisation_frame=PolarisationFrame('stokesI'))
def setUp(self): self.dir = './test_results' os.makedirs(self.dir, exist_ok=True) self.niter = 1000 self.lowcore = create_named_configuration('LOWBD2-CORE') self.nchan = 5 self.times = (numpy.pi / 12.0) * numpy.linspace(-3.0, 3.0, 7) self.frequency = numpy.linspace(0.9e8, 1.1e8, self.nchan) self.channel_bandwidth = numpy.array(self.nchan * [self.frequency[1] - self.frequency[0]]) self.phasecentre = SkyCoord(ra=+0.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000') self.vis = create_visibility(self.lowcore, self.times, self.frequency, self.channel_bandwidth, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame('stokesI')) self.vis.data['vis'] *= 0.0 # Create model self.test_model = create_low_test_image_from_gleam(npixel=512, cellsize=0.001, phasecentre=self.vis.phasecentre, frequency=self.frequency, channel_bandwidth=self.channel_bandwidth) beam = create_low_test_beam(self.test_model) export_image_to_fits(beam, "%s/test_deconvolve_msmfsclean_beam.fits" % self.dir) self.test_model.data *= beam.data export_image_to_fits(self.test_model, "%s/test_deconvolve_msmfsclean_model.fits" % self.dir) self.vis = predict_2d(self.vis, self.test_model) assert numpy.max(numpy.abs(self.vis.vis)) > 0.0 self.model = create_image_from_visibility(self.vis, npixel=512, cellsize=0.001, polarisation_frame=PolarisationFrame('stokesI')) self.dirty, sumwt = invert_2d(self.vis, self.model) self.psf, sumwt = invert_2d(self.vis, self.model, dopsf=True) export_image_to_fits(self.dirty, "%s/test_deconvolve_msmfsclean_dirty.fits" % self.dir) export_image_to_fits(self.psf, "%s/test_deconvolve_msmfsclean_psf.fits" % self.dir) window = numpy.ones(shape=self.model.shape, dtype=numpy.bool) window[..., 129:384, 129:384] = True self.innerquarter = create_image_from_array(window, self.model.wcs)
def setUp(self): self.lowcore = create_named_configuration('LOWBD2-CORE') self.times = numpy.linspace(-300.0, 300.0, 11) * numpy.pi / 43200.0 self.frequency = numpy.linspace(1e8, 1.5e9, 7) self.channel_bandwidth = numpy.array(7 * [self.frequency[1] - self.frequency[0]]) self.phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
def setUp(self): self.dir = './test_results' os.makedirs(self.dir, exist_ok=True) self.lowcore = create_named_configuration('LOWBD2-CORE') self.times = (numpy.pi / (12.0)) * numpy.linspace(-3.0, 3.0, 7) self.frequency = numpy.array([1e8]) self.channel_bandwidth = numpy.array([1e6]) self.phasecentre = SkyCoord(ra=+180.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000') self.vis = create_visibility(self.lowcore, self.times, self.frequency, channel_bandwidth=self.channel_bandwidth, phasecentre=self.phasecentre, weight=1.0, polarisation_frame=PolarisationFrame('stokesI')) self.vis.data['vis'] *= 0.0 # Create model self.test_model = create_test_image(cellsize=0.001, phasecentre=self.vis.phasecentre, frequency=self.frequency) self.vis = predict_2d(self.vis, self.test_model) assert numpy.max(numpy.abs(self.vis.vis)) > 0.0 self.model = create_image_from_visibility(self.vis, npixel=512, cellsize=0.001, polarisation_frame=PolarisationFrame('stokesI')) self.dirty, sumwt = invert_2d(self.vis, self.model) self.psf, sumwt = invert_2d(self.vis, self.model, dopsf=True) window = numpy.zeros(shape=self.model.shape, dtype=numpy.bool) window[..., 129:384, 129:384] = True self.innerquarter = create_image_from_array(window, self.model.wcs)
