Transform reference catalog sky positions (self.all_radec)
        to reference tangent plane (self.wcs) to create output X,Y positions.
861
def transformToRef(self):
        
        if 'refxyunits' in self.pars and self.pars['refxyunits'] == 'pixels':
            log.info('Creating RA/Dec positions for reference sources...')
            self.outxy = np.column_stack([self.all_radec[0][:,np.newaxis],self.all_radec[1][:,np.newaxis]])
            skypos = self.wcs.wcs_pix2world(self.all_radec[0],self.all_radec[1],self.origin)
            self.all_radec[0] = skypos[0]
            self.all_radec[1] = skypos[1]
        else:
            log.info('Converting RA/Dec positions of reference sources from "%s" to '%self.name+
                        'X,Y positions in reference WCS...')
            self.refWCS = self.wcs
            outxy = self.wcs.wcs_world2pix(self.all_radec[0],self.all_radec[1],self.origin)
            # convert outxy list to a Nx2 array
            self.outxy = np.column_stack([outxy[0][:,np.newaxis],outxy[1][:,np.newaxis]])
478726
def xy2rd(self,pos):
        
        if self.ctype1.find('TAN') < 0 or self.ctype2.find('TAN') < 0:
            print('XY2RD only supported for TAN projections.')
            raise TypeError

        if isinstance(pos,N.ndarray):
            # If we are working with an array of positions,
            # point to just X and Y values
            posx = pos[:,0]
            posy = pos[:,1]
        else:
            # Otherwise, we are working with a single X,Y tuple
            posx = pos[0]
            posy = pos[1]

        xi = self.cd11 * (posx - self.crpix1) + self.cd12 * (posy - self.crpix2)
        eta = self.cd21 * (posx - self.crpix1) + self.cd22 * (posy - self.crpix2)

        xi = DEGTORAD(xi)
        eta = DEGTORAD(eta)
        ra0 = DEGTORAD(self.crval1)
        dec0 = DEGTORAD(self.crval2)

        ra = N.arctan((xi / (N.cos(dec0)-eta*N.sin(dec0)))) + ra0
        dec = N.arctan( ((eta*N.cos(dec0)+N.sin(dec0)) /
                (N.sqrt((N.cos(dec0)-eta*N.sin(dec0))**2 + xi**2))) )

        ra = RADTODEG(ra)
        dec = RADTODEG(dec)
        ra = DIVMOD(ra, 360.)

        # Otherwise, just return the RA,Dec tuple.
        return ra,dec
1066579
def create_ellipse(self,xcen,ycen,a,b,ang,resolution=40.0):
        """"""

        import math
        e1=[]
        e2=[]
        ang=ang-math.radians(90)
        for i in range(0,int(resolution)+1):
            x=(-1*a+2*a*float(i)/resolution)
            y=1-(x/a)**2
            if y < 1E-6:
                y=1E-6
            y=math.sqrt(y)*b
            ptv=self.p2c((x*math.cos(ang)+y*math.sin(ang)+xcen,y*math.cos(ang)-x*math.sin(ang)+ycen))
            y=-1*y
            ntv=self.p2c((x*math.cos(ang)+y*math.sin(ang)+xcen,y*math.cos(ang)-x*math.sin(ang)+ycen))
            e1.append(ptv)
            e2.append(ntv)
        e2.reverse()
        e1.extend(e2)
        self.create_line(e1,fill='red',width=1)
613925
def surface(x, y, z):
    
    data = [go.Surface(x=x, y=y, z=z)]
    return Chart(data=data)
730783
def _make_wcsgeom_from_config(config):
        """"""

        binning = config['binning']
        binsz = binning['binsz']
        coordsys = binning.get('coordsys', 'GAL')
        roiwidth = binning['roiwidth']
        proj = binning.get('proj', 'AIT')
        ra = config['selection']['ra']
        dec = config['selection']['dec']
        npix = int(np.round(roiwidth / binsz))
        skydir = SkyCoord(ra * u.deg, dec * u.deg)

        wcsgeom = WcsGeom.create(npix=npix, binsz=binsz,
                                 proj=proj, coordsys=coordsys,
                                 skydir=skydir)
        return wcsgeom
511776
def raster_to_projection_coords(self, pixel_x, pixel_y):
        
        h_px_py = np.array([1, pixel_x, pixel_y])
        gt = np.array([[1, 0, 0], self.geotransform[0:3], self.geotransform[3:6]])
        arr = np.inner(gt, h_px_py)
        return arr[2], arr[1]
175276
def gprmc_to_degdec(lat, latDirn, lng, lngDirn):
    """"""
    x = float(lat[0:2]) + float(lat[2:]) / 60
    y = float(lng[0:3]) + float(lng[3:]) / 60

    if latDirn == 'S':
        x = -x
    if lngDirn == 'W':
        y = -y

    return x, y
1231843
function() {
            if (this.polygons.length === 0) {
                var m = new CSG.Matrix4x4(); // unity
                return [m,m];
            } else {
                // get a list of unique planes in the CSG:
                var csg = this.canonicalized();
                var planemap = {};
                csg.polygons.map(function(polygon) {
                    planemap[polygon.plane.getTag()] = polygon.plane;
                });
                // try each plane in the CSG and find the plane that, when we align it flat onto z=0,
                // gives the least height in z-direction.
                // If two planes give the same height, pick the plane that originally had a normal closest
                // to [0,0,-1].
                var xvector = new CSG.Vector3D(1, 0, 0);
                var yvector = new CSG.Vector3D(0, 1, 0);
                var zvector = new CSG.Vector3D(0, 0, 1);
                var z0connectorx = new CSG.Connector([0, 0, 0], [0, 0, -1], xvector);
                var z0connectory = new CSG.Connector([0, 0, 0], [0, 0, -1], yvector);
                var isfirst = true;
                var minheight = 0;
                var maxdotz = 0;
                var besttransformation, bestinversetransformation;
                for (var planetag in planemap) {
                    var plane = planemap[planetag];
                    var pointonplane = plane.normal.times(plane.w);
                    var transformation, inversetransformation;
                    // We need a normal vecrtor for the transformation
                    // determine which is more perpendicular to the plane normal: x or y?
                    // we will align this as much as possible to the x or y axis vector
                    var xorthogonality = plane.normal.cross(xvector).length();
                    var yorthogonality = plane.normal.cross(yvector).length();
                    if (xorthogonality > yorthogonality) {
                        // x is better:
                        var planeconnector = new CSG.Connector(pointonplane, plane.normal, xvector);
                        transformation = planeconnector.getTransformationTo(z0connectorx, false, 0);
                        inversetransformation = z0connectorx.getTransformationTo(planeconnector, false, 0);
                    } else {
                        // y is better:
                        var planeconnector = new CSG.Connector(pointonplane, plane.normal, yvector);
                        transformation = planeconnector.getTransformationTo(z0connectory, false, 0);
                        inversetransformation = z0connectory.getTransformationTo(planeconnector, false, 0);
                    }
                    var transformedcsg = csg.transform(transformation);
                    var dotz = -plane.normal.dot(zvector);
                    var bounds = transformedcsg.getBounds();
                    var zheight = bounds[1].z - bounds[0].z;
                    var isbetter = isfirst;
                    if (!isbetter) {
                        if (zheight < minheight) {
                            isbetter = true;
                        } else if (zheight == minheight) {
                            if (dotz > maxdotz) isbetter = true;
                        }
                    }
                    if (isbetter) {
                        // translate the transformation around the z-axis and onto the z plane:
                        var translation = new CSG.Vector3D([-0.5 * (bounds[1].x + bounds[0].x), -0.5 * (bounds[1].y + bounds[0].y), -bounds[0].z]);
                        transformation = transformation.multiply(CSG.Matrix4x4.translation(translation));
                        inversetransformation = CSG.Matrix4x4.translation(translation.negated()).multiply(inversetransformation);
                        minheight = zheight;
                        maxdotz = dotz;
                        besttransformation = transformation;
                        bestinversetransformation = inversetransformation;
                    }
                    isfirst = false;
                }
                return [besttransformation, bestinversetransformation];
            }
        }
561695
def transform_to(target_datum)
      return self if datum == target_datum
      if helmert = Osgb::Helmert[:"#{self.datum}_to_#{target_datum}"]
        cartesian_coordinates = Osgb::Ellipsoid[self.datum].polar_to_cartesian(@lat.to_radians,@lng.to_radians)
        transformed = helmert.transform(*cartesian_coordinates)
        phi, lambda = Osgb::Ellipsoid[target_datum].cartesian_to_polar(*transformed)
        self.class.new(phi.to_degrees, lambda.to_degrees, target_datum, precision)
      else
        raise Osgb::TransformationError, "Missing helmert transformation for #{self.datum} to #{target_datum}"
      end
    end
711225
def _add_projection_coords(ds, prod_desc, proj_var, dx, dy):
    """"""
    proj = cf_to_proj(proj_var)

    # Get projected location of lower left point
    x0, y0 = proj(prod_desc.lo1, prod_desc.la1)

    # Coordinate variable for x
    ds.createDimension('x', prod_desc.nx)
    x_var = ds.createVariable('x', np.float64, dimensions=('x',))
    x_var.units = 'm'
    x_var.long_name = 'x coordinate of projection'
    x_var.standard_name = 'projection_x_coordinate'
    x_var[:] = x0 + np.arange(prod_desc.nx) * (1000. * dx)

    # Now y
    ds.createDimension('y', prod_desc.ny)
    y_var = ds.createVariable('y', np.float64, dimensions=('y',))
    y_var.units = 'm'
    y_var.long_name = 'y coordinate of projection'
    y_var.standard_name = 'projection_y_coordinate'

    # Need to flip y because we calculated from the lower left corner, but the raster data
    # is stored with top row first.
    y_var[::-1] = y0 + np.arange(prod_desc.ny) * (1000. * dy)

    # Get the two-D lon,lat grid as well
    x, y = np.meshgrid(x_var[:], y_var[:])
    lon, lat = proj(x, y, inverse=True)
    lon_var = ds.createVariable('lon', np.float64, dimensions=('y', 'x'), wrap_array=lon)
    lon_var.long_name = 'longitude'
    lon_var.units = 'degrees_east'

    lat_var = ds.createVariable('lat', np.float64, dimensions=('y', 'x'), wrap_array=lat)
    lat_var.long_name = 'latitude'
    lat_var.units = 'degrees_north'

    ds.img_extent = (x_var[:].min(), x_var[:].max(), y_var[:].min(), y_var[:].max())
630120
def reference_source_point(self):
        
        xref = isinstance(self.xref, Quantity) and self.xref.value or self.xref
        yref = isinstance(self.yref, Quantity) and self.yref.value or self.yref

        return xref + self.x_ref_offset, yref + self.y_ref_offset
1121223
def _map_from_binaries(self, eopatch, dst_shape, request_data):
        
        if self.feature_name in eopatch[self.feature_type]:
            raster = eopatch[self.feature_type][self.feature_name].squeeze()
        else:
            raster = np.ones(dst_shape, dtype=self.raster_dtype) * self.no_data_val

        new_raster = self._reproject(eopatch, self._to_binary_mask(request_data))

        # update raster
        raster[new_raster != 0] = new_raster[new_raster != 0]

        return raster
384502
def to_sky(self, wcs, origin=_DEFAULT_WCS_ORIGIN, mode=_DEFAULT_WCS_MODE):
        
        return SkyCoord.from_pixel(
            xp=self.x, yp=self.y, wcs=wcs,
            origin=origin, mode=mode,
        )
39225
def makeHDU(self, mag_1, mag_err_1, mag_2, mag_err_2, lon, lat, mc_source_id):
        
        if self.config['catalog']['coordsys'].lower() == 'cel' \
           and self.config['coords']['coordsys'].lower() == 'gal':
            lon, lat = ugali.utils.projector.gal2cel(lon, lat)
        elif self.config['catalog']['coordsys'].lower() == 'gal' \
           and self.config['coords']['coordsys'].lower() == 'cel':
            lon, lat = ugali.utils.projector.cel2gal(lon, lat)

        columns = [
            pyfits.Column(name=self.config['catalog']['objid_field'],
                          format = 'D',array = np.arange(len(lon))),
            pyfits.Column(name=self.config['catalog']['lon_field'],
                          format = 'D',array = lon),
            pyfits.Column(name = self.config['catalog']['lat_field'],          
                          format = 'D',array = lat), 
            pyfits.Column(name = self.config['catalog']['mag_1_field'],        
                          format = 'E',array = mag_1),
            pyfits.Column(name = self.config['catalog']['mag_err_1_field'],    
                          format = 'E',array = mag_err_1),
            pyfits.Column(name = self.config['catalog']['mag_2_field'],        
                          format = 'E',array = mag_2),
            pyfits.Column(name = self.config['catalog']['mag_err_2_field'],    
                          format = 'E',array = mag_err_2),
            pyfits.Column(name = self.config['catalog']['mc_source_id_field'], 
                          format = 'I',array = mc_source_id),
        ]

        hdu = pyfits.new_table(columns)
        return hdu
816222
def do_driz(insci, input_wcs, inwht,
            output_wcs, outsci, outwht, outcon,
            expin, in_units, wt_scl,
            wcslin_pscale=1.0,uniqid=1, pixfrac=1.0, kernel='square',
            fillval="INDEF", stepsize=10,wcsmap=None):
    
    # Insure that the fillval parameter gets properly interpreted for use with tdriz
    if util.is_blank(fillval):
        fillval = 'INDEF'
    else:
        fillval = str(fillval)

    if in_units == 'cps':
        expscale = 1.0
    else:
        expscale = expin

    # Compute what plane of the context image this input would
    # correspond to:
    planeid = int((uniqid-1) / 32)

    # Check if the context image has this many planes
    if outcon.ndim == 3:
        nplanes = outcon.shape[0]
    elif outcon.ndim == 2:
        nplanes = 1
    else:
        nplanes = 0

    if nplanes <= planeid:
        raise IndexError("Not enough planes in drizzle context image")

    # Alias context image to the requested plane if 3d
    if outcon.ndim == 2:
        outctx = outcon
    else:
        outctx = outcon[planeid]

    pix_ratio = output_wcs.pscale/wcslin_pscale

    if wcsmap is None and cdriz is not None:
        log.info('Using WCSLIB-based coordinate transformation...')
        log.info('stepsize = %s' % stepsize)
        mapping = cdriz.DefaultWCSMapping(
            input_wcs, output_wcs,
            input_wcs.pixel_shape[0], input_wcs.pixel_shape[1],
            stepsize
        )
    else:
        #
        ##Using the Python class for the WCS-based transformation
        #
        # Use user provided mapping function
        log.info('Using coordinate transformation defined by user...')
        if wcsmap is None:
            wcsmap = wcs_functions.WCSMap
        wmap = wcsmap(input_wcs,output_wcs)
        mapping = wmap.forward

    _shift_fr = 'output'
    _shift_un = 'output'
    ystart = 0
    nmiss = 0
    nskip = 0
    #
    # This call to 'cdriz.tdriz' uses the new C syntax
    #
    _dny = insci.shape[0]
    # Call 'drizzle' to perform image combination
    if insci.dtype > np.float32:
        #WARNING: Input array recast as a float32 array
        insci = insci.astype(np.float32)

    _vers,nmiss,nskip = cdriz.tdriz(insci, inwht, outsci, outwht,
        outctx, uniqid, ystart, 1, 1, _dny,
        pix_ratio, 1.0, 1.0, 'center', pixfrac,
        kernel, in_units, expscale, wt_scl,
        fillval, nmiss, nskip, 1, mapping)

    if nmiss > 0:
        log.warning('! %s points were outside the output image.' % nmiss)
    if nskip > 0:
        log.debug('! Note, %s input lines were skipped completely.' % nskip)

    return _vers
766013
def make_fermipy_roi_model_from_catalogs(cataloglist):
        
        data = dict(catalogs=cataloglist,
                    src_roiwidth=360.)
        return roi_model.ROIModel(data, skydir=SkyCoord(0.0, 0.0, unit='deg'))
1259465
def rd2xy(input,ra=None,dec=None,coordfile=None,colnames=None,
            precision=6,output=None,verbose=True):
    
    single_coord = False
    if coordfile is not None:
        if colnames in blank_list:
            colnames = ['c1','c2']
        elif isinstance(colnames,type('a')):
            colnames = colnames.split(',')
        # convert input file coordinates to lists of decimal degrees values
        xlist,ylist = tweakutils.readcols(coordfile,cols=colnames)
    else:
        if isinstance(ra,np.ndarray):
            ralist = ra.tolist()
            declist = dec.tolist()
        elif not isinstance(ra, list):
            ralist = [ra]
            declist = [dec]
        else:
            ralist = ra
            declist = dec
        xlist  = [0]*len(ralist)
        ylist = [0]*len(ralist)
        if len(xlist) == 1:
            single_coord = True
        for i,(r,d) in enumerate(zip(ralist,declist)):
            # convert input value into decimal degrees value
            xval,yval = tweakutils.parse_skypos(r,d)
            xlist[i] = xval
            ylist[i] = yval

    # start by reading in WCS+distortion info for input image
    inwcs = wcsutil.HSTWCS(input)
    if inwcs.wcs.is_unity():
        print("####\nNo valid WCS found in {}.\n  Results may be invalid.\n####\n".format(input))

    # Now, convert pixel coordinates into sky coordinates
    try:
        outx,outy = inwcs.all_world2pix(xlist,ylist,1)
    except RuntimeError:
        outx,outy = inwcs.wcs_world2pix(xlist,ylist,1)

    # add formatting based on precision here...
    xstr = []
    ystr = []
    fmt = "%."+repr(precision)+"f"
    for x,y in zip(outx,outy):
        xstr.append(fmt%x)
        ystr.append(fmt%y)

    if verbose or (not verbose and util.is_blank(output)):
        print ('# Coordinate transformations for ',input)
        print('# X      Y         RA             Dec\n')
        for x,y,r,d in zip(xstr,ystr,xlist,ylist):
            print("%s  %s    %s  %s"%(x,y,r,d))

    # Create output file, if specified
    if output:
        f = open(output,mode='w')
        f.write("# Coordinates converted from %s\n"%input)
        for x,y in zip(xstr,ystr):
            f.write('%s    %s\n'%(x,y))
        f.close()
        print('Wrote out results to: ',output)

    if single_coord:
        outx = outx[0]
        outy = outy[0]

    return outx, outy