我们从Python开源项目中，提取了以下25个代码示例，用于说明如何使用astropy.units.rad()。
def calc_rate(redshift): """Calculate the rate, kpc/px.""" # Init # Hubble constant at z = 0 H0 = 71.0 # Omega0, total matter density Om0 = 0.27 # Cosmo cosmo = FlatLambdaCDM(H0=H0, Om0=Om0) # Angular diameter distance, [Mpc] dA = cosmo.angular_diameter_distance(redshift) # Resolution of Chandra rate_px = 0.492 * au.arcsec # [arcmin/px] rate_px_rad = rate_px.to(au.rad) # rate rate = dA.value * rate_px_rad.value # [Mpc/px] return rate
def test_phase_rotation(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) # Predict visibilities with new phase centre independently ha_diff = -(self.compabsdirection.ra - self.phasecentre.ra).to(u.rad).value vispred = create_visibility(self.lowcore, self.times + ha_diff, self.frequency, channel_bandwidth=self.channel_bandwidth, phasecentre=self.compabsdirection, weight=1.0, polarisation_frame=PolarisationFrame("stokesIQUV")) vismodel2 = predict_skycomponent_visibility(vispred, self.comp) # Should yield the same results as rotation rotatedvis = phaserotate_visibility(self.vismodel, newphasecentre=self.compabsdirection, tangent=False) assert_allclose(rotatedvis.vis, vismodel2.vis, rtol=1e-7) assert_allclose(rotatedvis.uvw, vismodel2.uvw, rtol=1e-7)
def NES_healpixels(nside=set_default_nside(), width=15, dec_min=0., fill_gap=True): """ Define the North Ecliptic Spur region. Return a healpix map with NES pixels as 1. """ ra, dec = ra_dec_hp_map(nside=nside) result = np.zeros(ra.size) coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad) eclip_lat = coord.barycentrictrueecliptic.lat.radian good = np.where((np.abs(eclip_lat) <= np.radians(width)) & (dec > dec_min)) result[good] += 1 if fill_gap: good = np.where((dec > np.radians(dec_min)) & (ra < np.radians(180)) & (dec < np.radians(width))) result[good] = 1 return result
def galactic_plane_healpixels(nside=set_default_nside(), center_width=10., end_width=4., gal_long1=70., gal_long2=290.): """ Define the Galactic Plane region. Return a healpix map with GP pixels as 1. """ ra, dec = ra_dec_hp_map(nside=nside) result = np.zeros(ra.size) coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad) g_long, g_lat = coord.galactic.l.radian, coord.galactic.b.radian good = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) < np.radians(center_width))) result[good] += 1 good = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) < np.radians(center_width))) result[good] += 1 # Add tapers slope = -(np.radians(center_width)-np.radians(end_width))/(np.radians(gal_long1)) lat_limit = slope*g_long+np.radians(center_width) outside = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) > np.abs(lat_limit))) result[outside] = 0 slope = (np.radians(center_width)-np.radians(end_width))/(np.radians(360. - gal_long2)) b = np.radians(center_width)-np.radians(360.)*slope lat_limit = slope*g_long+b outside = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) > np.abs(lat_limit))) result[outside] = 0 return result
def is_NES(nside=64, width=15, dec_min=0., fill_gap=True): """ Define the North Ecliptic Spur region. Return a healpix map with NES pixels as true. """ ra, dec = ra_dec_hp_map(nside=nside) result = np.zeros(ra.size) coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad) eclip_lat = coord.barycentrictrueecliptic.lat.radian good = np.where((np.abs(eclip_lat) <= np.radians(width)) & (dec > dec_min)) result[good] += 1 if fill_gap: good = np.where((dec > np.radians(dec_min)) & (ra < np.radians(180)) & (dec < np.radians(width))) result[good] = 1 NES_indx = (result==1) return NES_indx
def rho_as_angle(rho, eph): """Projected linear distance to angular distance. Parameters ---------- rho : `~astropy.units.Quantity` Projected distance in units of length. eph : dictionary-like or `~sbpy.data.Ephem` Ephemerides; requires geocentric distance as `delta`. Returns ------- rho_l : `~astropy.units.Quantity` """ if rho.unit.is_equivalent(u.m): rho_a = np.arctan(rho / eph['delta'].to(u.m)) else: assert rho.unit.is_equivalent(u.rad) rho_a = rho return rho_a
def rho_as_length(rho, eph): """Angular distance to projected linear distance. Parameters ---------- rho : `~astropy.units.Quantity` Projected distance in units of angle. eph : dictionary-like or `~sbpy.data.Ephem` Ephemerides; requires geocentric distance as `delta`. Returns ------- rho_l : `~astropy.units.Quantity` """ if rho.unit.is_equivalent(u.rad): rho_l = eph['delta'].to(u.m) * np.tan(rho) else: assert rho.unit.is_equivalent(u.m) rho_l = rho return rho_l
def create_configuration_from_file(antfile: str, name: str = None, location: EarthLocation = None, mount: str = 'altaz', names: str = "%d", frame: str = 'local', diameter=35.0, meta: dict = None, rmax=None, **kwargs) -> Configuration: """ Define from a file :param names: :param antfile: Antenna file name :param name: Name of array e.g. 'LOWBD2' :param location: :param mount: mount type: 'altaz', 'xy' :param frame: 'local' | 'global' :param diameter: Effective diameter of station or antenna :param meta: Any meta info :return: Configuration """ antxyz = numpy.genfromtxt(antfile, delimiter=",") assert antxyz.shape[1] == 3, ("Antenna array has wrong shape %s" % antxyz.shape) if frame == 'local': latitude = location.geodetic[1].to(u.rad).value antxyz = xyz_at_latitude(antxyz, latitude) if rmax is not None: lantxyz = antxyz - numpy.average(antxyz, axis=0) r = numpy.sqrt(lantxyz[:, 0] ** 2 + lantxyz[:, 1] ** 2 + lantxyz[:, 2] ** 2) antxyz = antxyz[r < rmax] log.debug('create_configuration_from_file: Maximum radius %.1f m includes %d antennas/stations' % (rmax, antxyz.shape[0])) nants = antxyz.shape[0] anames = [names % ant for ant in range(nants)] mounts = numpy.repeat(mount, nants) fc = Configuration(location=location, names=anames, mount=mounts, xyz=antxyz, frame=frame, diameter=diameter) return fc
def kpc_per_arcsec(self, z): """ Calculate the transversal length (unit: kpc) corresponding to 1 arcsec at the *angular diameter distance* of z. """ dist_kpc = self.angular_diameter_distance(z, unit="kpc") return dist_kpc * au.arcsec.to(au.rad)
def kpc_per_pix(self, z): """ Calculate the transversal length (unit: kpc) corresponding to 1 ACIS pixel (i.e., 0.492 arcsec) at the *angular diameter distance* of z. """ pix = ACIS.pixel2arcsec * au.arcsec dist_kpc = self.angular_diameter_distance(z, unit="kpc") return dist_kpc * pix.to(au.rad).value
def find_contours(data, segments, sigma_level): """ This function ... :param data: :param segments: :param sigma_level: :return: """ # Initialize a list for the contours contours = [] # Get the segment properties # Since there is only one segment in the source.mask (the center segment), the props # list contains only one entry (one galaxy) properties_list = source_properties(data, segments) for properties in properties_list: # Obtain the position, orientation and extent position = Position(properties.xcentroid.value, properties.ycentroid.value) a = properties.semimajor_axis_sigma.value * sigma_level b = properties.semiminor_axis_sigma.value * sigma_level angle = properties.orientation.value # in radians angle = Angle(angle, u.rad) radius = Extent(a, b) meta = {"text": str(properties.label)} # Create the contour contours.append(Ellipse(position, radius, angle, meta=meta)) # Return the contours return contours # -----------------------------------------------------------------
def find_contour(box, mask, sigma_level): """ This function ... :param box: :param mask: :param sigma_level: :return: """ props = source_properties(box, mask) #tbl = properties_table(props) x_shift = box.x_min y_shift = box.y_min # Since there is only one segment in the self.source.mask (the center segment), the props # list contains only one entry (one galaxy) if len(props) == 0: return None properties = props[0] # Obtain the position, orientation and extent position = Position(properties.xcentroid.value + x_shift, properties.ycentroid.value + y_shift) a = properties.semimajor_axis_sigma.value * sigma_level b = properties.semiminor_axis_sigma.value * sigma_level angle = properties.orientation.value # in radians angle = Angle(angle, u.rad) radius = Extent(a, b) # Create and return the elliptical contour return Ellipse(position, radius, angle) # -----------------------------------------------------------------
def calculate(self): ephem_location = ephem.Observer() ephem_location.lat = self.location.latitude.to(u.rad) / u.rad ephem_location.lon = self.location.longitude.to(u.rad) / u.rad ephem_location.elevation = self.location.height / u.meter ephem_location.date = ephem.Date(self.time.datetime) if self.data is None: self.alt = Latitude([], unit=u.deg) self.az = Longitude([], unit=u.deg) self.names = Column([], dtype=np.str) self.vmag = Column([]) else: ra = Longitude((self.data['RAh'], self.data['RAm'], self.data['RAs']), u.h) dec = Latitude((np.core.defchararray.add(self.data['DE-'], self.data['DEd'].astype(str)).astype(int), self.data['DEm'], self.data['DEs']), u.deg) c = SkyCoord(ra, dec, frame='icrs') altaz = c.transform_to(AltAz(obstime=self.time, location=self.location)) self.alt = altaz.alt self.az = altaz.az self.names = self.data['Name'] self.vmag = self.data['Vmag'] for ephemeris in self.ephemerides: ephemeris.compute(ephem_location) self.vmag = self.vmag.insert(0, ephemeris.mag) self.alt = self.alt.insert(0, (ephemeris.alt.znorm * u.rad).to(u.deg)) self.az = self.az.insert(0, (ephemeris.az * u.rad).to(u.deg)) self.names = self.names.insert(0, ephemeris.name) return self.names, self.vmag, self.alt, self.az
def selectSrc(self): """ Filter a given source, running gtselect """ # Do we have to deal with a FITS file or an ASCII list of FITS file ? allskyext = os.path.splitext(self.allsky)[1] if allskyext in [".fit", ".fits"]: fermi.filter['infile'] = self.allsky else: fermi.filter['infile'] = '@%s' % self.allsky if self.daily: outfile = self.workDir + '/' + str(self.src) + '_daily.fits' else: outfile = self.workDir + '/' + str(self.src) + '.fits' fermi.filter['outfile'] = outfile # If outfile already exists, we don't do anything if os.path.isfile(outfile): return True fermi.filter['ra'] = self.ra fermi.filter['dec'] = self.dec fermi.filter['rad'] = self.roi fermi.filter['emin'] = self.emin fermi.filter['emax'] = self.emax fermi.filter['tmin'] = self.tstart fermi.filter['tmax'] = self.tstop fermi.filter['zmax'] = self.zmax fermi.filter['evclass'] = 128 logging.info('Running gtselect') fermi.filter.run()
def exposure(self, gamma=None): """ Compute exposure on source src, to add a flux column for the photometric light curve. Warning: the input file is modified in place, with an additional exposure column added to the file ! """ if self.daily: infile = self.workDir + '/' + str(self.src) + '_daily_lc.fits' srcmdl = self.workDir + '/' + str(self.src) + '_daily.xml' else: infile = self.workDir + '/' + str(self.src) + '_lc.fits' srcmdl = self.workDir + '/' + str(self.src) + '.xml' # If infile already contains an EXPOSURE column, we don't do anything hdu = fits.open(infile) if hdu[1].header.get('TTYPE5') == 'EXPOSURE': return True scfile = self.spacecraft irfs = 'P8R2_SOURCE_V6' rad = str(self.roi) options = 'infile=' + infile + ' scfile=' + scfile + ' irfs=' + irfs + ' rad=' + rad if self.fglName is not None: target = self.fglName.replace('3FGLJ', '3FGL J') logging.debug('exposure: target=%s', target) options += ' srcmdl=' + srcmdl + ' target="' + target + '"' else: options += ' srcmdl="none" specin=' + str(gamma) cmd = 'time -p ' + self.fermiDir + '/bin/gtexposure ' + options logging.info('Running gtexposure') os.system(cmd)
def stupidFast_altAz2RaDec(alt, az, lat, lon, mjd): """ Convert alt, az to RA, Dec without taking into account abberation, precesion, diffraction, ect. Parameters ---------- alt : numpy.array Altitude, same length as `ra` and `dec`. Radians. az : numpy.array Azimuth, same length as `ra` and `dec`. Must be same length as `alt`. Radians. lat : float Latitude of the observatory in radians. lon : float Longitude of the observatory in radians. mjd : float Modified Julian Date. Returns ------- ra : array_like RA, in radians. dec : array_like Dec, in radians. """ lmst, last = calcLmstLast(mjd, lon) lmst = lmst/12.*np.pi # convert to rad sindec = np.sin(lat)*np.sin(alt) + np.cos(lat)*np.cos(alt)*np.cos(az) sindec = inrange(sindec) dec = np.arcsin(sindec) ha = np.arctan2(-np.sin(az)*np.cos(alt), -np.cos(az)*np.sin(lat)*np.cos(alt)+np.sin(alt)*np.cos(lat)) ra = (lmst-ha) raneg = np.where(ra < 0) ra[raneg] = ra[raneg] + 2.*np.pi raover = np.where(ra > 2.*np.pi) ra[raover] -= 2.*np.pi return ra, dec
def is_GP(nside=64, center_width=10., end_width=4., gal_long1=70., gal_long2=290.): """ Define the Galactic Plane region. Return a healpix map with GP pixels as true. """ ra, dec = ra_dec_hp_map(nside=nside) result = np.zeros(ra.size) coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad) g_long, g_lat = coord.galactic.l.radian, coord.galactic.b.radian good = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) < np.radians(center_width))) result[good] += 1 good = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) < np.radians(center_width))) result[good] += 1 # Add tapers slope = -(np.radians(center_width)-np.radians(end_width))/(np.radians(gal_long1)) lat_limit = slope*g_long+np.radians(center_width) outside = np.where((g_long < np.radians(gal_long1)) & (np.abs(g_lat) > np.abs(lat_limit))) result[outside] = 0 slope = (np.radians(center_width)-np.radians(end_width))/(np.radians(360. - gal_long2)) b = np.radians(center_width)-np.radians(360.)*slope lat_limit = slope*g_long+b outside = np.where((g_long > np.radians(gal_long2)) & (np.abs(g_lat) > np.abs(lat_limit))) result[outside] = 0 GP_inx =(result==1) return GP_inx
def xarray(shape, yx=[0, 0], th=None): """2D array of x values. Parameters ---------- shape : tuple of int The shape of the resulting array, (y, x). yx : tuple of float Offset the array to align with this y, x center. th : Quanitity, optional Place the x-axis along this position angle, measured counterclockwise from the original x-axis. Returns ------- x : ndarray An array of x values. Examples -------- >>> from sbpy.imageanalysis.utils import xarray >>> x = xarray((10, 10)) >>> x[0, 3] 3 """ import numpy as np import astropy.units as u y, x = np.indices(shape)[-2:] y = y - yx[0] x = x - yx[1] if th is not None: x = x * np.cos(th.to(u.rad).value) + y * np.sin(th.to(u.rad).value) return x
def yarray(shape, yx=[0, 0], th=None): """2D array of y values. Parameters ---------- shape : tuple of int The shape of the resulting array, (y, x). yx : tuple of float Offset the array to align with this y, x center. th : Quantity, optional Place the y-axis along this position angle, measured counterclockwise from the original y-axis. Returns ------- y : ndarray An array of y values. >>> from sbpy.imageanalysis.utils import yarray >>> y = yarray((10, 10)) >>> y[3, 0] 3 """ import numpy as np import astropy.units as u y, x = np.indices(shape)[-2:] y = y - yx[0] x = x - yx[1] if th is not None: y = -x * np.sin(th.to(u.rad).value) + y * np.cos(th.to(u.rad).value) return y
def phaserotate_visibility(vis: Visibility, newphasecentre: SkyCoord, tangent=True, inverse=False) -> Visibility: """ Phase rotate from the current phase centre to a new phase centre If tangent is False the uvw are recomputed and the visibility phasecentre is updated. Otherwise only the visibility phases are adjusted :param vis: Visibility to be rotated :param newphasecentre: :param tangent: Stay on the same tangent plane? (True) :param inverse: Actually do the opposite :return: Visibility """ assert isinstance(vis, Visibility), "vis is not a Visibility: %r" % vis l, m, n = skycoord_to_lmn(newphasecentre, vis.phasecentre) # No significant change? if numpy.abs(n) > 1e-15: # Make a new copy newvis = copy_visibility(vis) phasor = simulate_point(newvis.uvw, l, m) nvis, npol = vis.vis.shape # TODO: Speed up (broadcast rules not obvious to me) if inverse: for i in range(nvis): for pol in range(npol): newvis.data['vis'][i, pol] *= phasor[i] else: for i in range(nvis): for pol in range(npol): newvis.data['vis'][i, pol] *= numpy.conj(phasor[i]) # To rotate UVW, rotate into the global XYZ coordinate system and back. We have the option of # staying on the tangent plane or not. If we stay on the tangent then the raster will # join smoothly at the edges. If we change the tangent then we will have to reproject to get # the results on the same image, in which case overlaps or gaps are difficult to deal with. if not tangent: if inverse: xyz = uvw_to_xyz(vis.data['uvw'], ha=-newvis.phasecentre.ra.rad, dec=newvis.phasecentre.dec.rad) newvis.data['uvw'][...] = \ xyz_to_uvw(xyz, ha=-newphasecentre.ra.rad, dec=newphasecentre.dec.rad)[...] else: # This is the original (non-inverse) code xyz = uvw_to_xyz(newvis.data['uvw'], ha=-newvis.phasecentre.ra.rad, dec=newvis.phasecentre.dec.rad) newvis.data['uvw'][...] = xyz_to_uvw(xyz, ha=-newphasecentre.ra.rad, dec=newphasecentre.dec.rad)[ ...] newvis.phasecentre = newphasecentre return newvis else: return vis
def mergeGTIfiles(self): """ Merge multiple GTI files when mergelongterm is True. Use gtselect. Assume the current workDir is longTerm/merged. """ # Create list of GTI files if not self.daily: listname = self.workDir + '/' + self.src + '_gti.list' else: listname = self.workDir + '/' + self.src + '_daily_gti.list' filelist = open(listname, 'w') list = [] if not self.daily: for file in glob.glob(self.workDir + '/../20????/' + self.src + '_gti.fits'): list.append(file) else: for file in glob.glob(self.workDir + '/../20????/' + self.src + '_daily_gti.fits'): list.append(file) # Sort the list of GTI files list = sorted(list) for item in list: filelist.write(item + '\n') filelist.close() fermi.filter['infile'] = '@' + listname if not self.daily: outfile = self.workDir + '/' + str(self.src) + '_gti.fits' else: outfile = self.workDir + '/' + str(self.src) + '_daily_gti.fits' fermi.filter['outfile'] = outfile # If outfile already exists, we re-create it if os.path.isfile(outfile): os.remove(outfile) fermi.filter['ra'] = self.ra fermi.filter['dec'] = self.dec fermi.filter['rad'] = self.roi fermi.filter['emin'] = self.emin fermi.filter['emax'] = self.emax fermi.filter['tmin'] = self.tstart fermi.filter['tmax'] = self.tstop fermi.filter['zmax'] = self.zmax fermi.filter['evclass'] = 128 logging.info('Running gtselect') fermi.filter.run()
def empty_observation(): """ Return a numpy array that could be a handy observation record XXX: Should this really be "empty visit"? Should we have "visits" made up of multple "observations" to support multi-exposure time visits? XXX-Could add a bool flag for "observed". Then easy to track all proposed observations. Could also add an mjd_min, mjd_max for when an observation should be observed. That way we could drop things into the queue for DD fields. XXX--might be nice to add a generic "sched_note" str field, to record any metadata that would be useful to the scheduler once it's observed. and/or observationID. Returns ------- numpy array Notes ----- The numpy fields have the following structure RA : float The Right Acension of the observation (center of the field) (Radians) dec : float Declination of the observation (Radians) mjd : float Modified Julian Date at the start of the observation (time shutter opens) exptime : float Total exposure time of the visit (seconds) filter : str The filter used. Should be one of u, g, r, i, z, y. rotSkyPos : float The rotation angle of the camera relative to the sky E of N (Radians) nexp : int Number of exposures in the visit. airmass : float Airmass at the center of the field FWHMeff : float The effective seeing FWHM at the center of the field. (arcsec) skybrightness : float The surface brightness of the sky background at the center of the field. (mag/sq arcsec) night : int The night number of the observation (days) """ names = ['RA', 'dec', 'mjd', 'exptime', 'filter', 'rotSkyPos', 'nexp', 'airmass', 'FWHMeff', 'FWHM_geometric', 'skybrightness', 'night', 'slewtime', 'fivesigmadepth', 'alt', 'az', 'clouds', 'moonAlt', 'sunAlt', 'note'] # units of rad, rad, days, seconds, string, radians (E of N?) types = [float, float, float, float, '|U1', float, int, float, float, float, float, int, float, float, float, float, float, float, float, '|U40'] result = np.zeros(1, dtype=list(zip(names, types))) return result
def stupidFast_RaDec2AltAz(ra, dec, lat, lon, mjd, lmst=None): """ Convert Ra,Dec to Altitude and Azimuth. Coordinate transformation is killing performance. Just use simple equations to speed it up and ignore abberation, precesion, nutation, nutrition, etc. Parameters ---------- ra : array_like RA, in radians. dec : array_like Dec, in radians. Must be same length as `ra`. lat : float Latitude of the observatory in radians. lon : float Longitude of the observatory in radians. mjd : float Modified Julian Date. Returns ------- alt : numpy.array Altitude, same length as `ra` and `dec`. Radians. az : numpy.array Azimuth, same length as `ra` and `dec`. Radians. """ if lmst is None: lmst, last = calcLmstLast(mjd, lon) lmst = lmst/12.*np.pi # convert to rad ha = lmst-ra sindec = np.sin(dec) sinlat = np.sin(lat) coslat = np.cos(lat) sinalt = sindec*sinlat+np.cos(dec)*coslat*np.cos(ha) sinalt = inrange(sinalt) alt = np.arcsin(sinalt) cosaz = (sindec-np.sin(alt)*sinlat)/(np.cos(alt)*coslat) cosaz = inrange(cosaz) az = np.arccos(cosaz) signflip = np.where(np.sin(ha) > 0) az[signflip] = 2.*np.pi-az[signflip] return alt, az
