How to use the tractor.sdss function in tractor

print rcf
        try:
            TI.extend([st.get_tractor_image(rcf[0], rcf[1], rcf[2], bandname,useMags=True) for bandname in bands])
            sources.append(st.get_tractor_sources(rcf[0], rcf[1], rcf[2],bandname,bands=bands))
        except:
            print "Sources not in data! %s" % (rcf)

    timg,info = TI[0]
    photocal = timg.getPhotoCal()

    wcs = timg.getWcs()
    lvl = logging.DEBUG
    logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)

    tims = [timg for timg,tinf in TI]
    tractor = st.SDSSTractor(tims)
    for source in sources:
        tractor.addSources(source)

    zr = np.array([-5.,+5.]) * info['skysig']

    print bands

    prefix = 'ngc4258'
    saveAll('initial-'+prefix, tractor,zr,flipBands,debug=True)
    plotInvvar('initial-'+prefix,tractor)
    assert(False)
    bright = None
    lowbright = 1000
    for timg,sources in zip(tims,sources):
        wcs = timg.getWcs()
        xtr,ytr = wcs.positionToPixel(RaDecPos(ra,dec))

def fp():
	band = 'i'
	ps = PlotSequence('fp')
	
	run, camcol, field = 5115, 5, 151
	bands = [band]
	roi = (1048,2048, 0,1000)
	tim,tinf = st.get_tractor_image_dr9(run, camcol, field, band,
										roi=roi, nanomaggies=True)
	ima = dict(interpolation='nearest', origin='lower',
			   extent=roi)
	zr2 = tinf['sky'] + tinf['skysig'] * np.array([-3, 100])
	#imb = ima.copy()
	#imb.update(vmin=tim.zr[0], vmax=tim.zr[1])
	imc = ima.copy()
	imc.update(norm=ArcsinhNormalize(mean=tinf['sky'], 
									 std=tinf['skysig']),
				vmin=zr2[0], vmax=zr2[1])


	

	# Match spectra with Abell catalog
	T = fits_table('a1656-spectro.fits', column_map={'class':'clazz'})

wcs = ba.FitsWcs(tan)

    data=np.zeros((imagew,imageh))
    invvar=np.ones((imagew,imageh))

    psf = ba.GaussianMixturePSF(1.,[0.,0.],np.array(1.)) #amp,mean,var
    
    skyobj = ba.ConstantSky(0.)
    zr = np.array([-5.,+5.])


    tims = []
    bands = ['u','g','r','i','z']
    for bandname in bands:
        photocal = st.SdssNanomaggiesPhotoCal(bandname)
        image = en.Image(data=data,invvar=invvar,sky=skyobj,psf=psf,wcs=wcs,photocal=photocal,name="Half-light %s" %bandname,zr=zr)
        tims.append(image)
    tractor = st.SDSSTractor(tims)
    tractor.addSources([CG])

    yg, xg = np.meshgrid(np.arange(imageh) - crpix2, np.arange(imagew) - crpix1)
    r2g = xg ** 2 + yg ** 2
    rlist_pix = np.exp(np.linspace(0.,np.log(0.5*imageh),64))
    rlist_arcsec = rlist_pix * pixscale * 3600.
    mimgs = tractor.getModelImages()
    r50s = []
    r90s = []
    concs = []
    expr50 = CG.shapeExp.re * np.sqrt(CG.shapeExp.ab)
    devr50 = CG.shapeDev.re * np.sqrt(CG.shapeDev.ab)

sumbright = sum([brightE.getMag(bandname)+brightD.getMag(bandname) for bandname in bands])
                if sumbright < lowbright:
                    print("GREATER")
                    lowBrightE = brightE
                    lowBrightD = brightD
                    lowShapeE = src.shapeExp
                    lowShapeD = src.shapeDev
            print "Removed:", src
            print xs,ys
            tractor.removeSource(src)

#    saveBands('removed-'+prefix, tractor,zr,flipBands,debug=True)
    plotInvvar('removed-'+prefix,tractor)


    CG = st.CompositeGalaxy(RaDecPos(ra,dec),lowBrightE,lowShapeE,lowBrightD,lowShapeD)
    print CG
    tractor.addSource(CG)


 #   saveBands('added-'+prefix,tractor,zr,flipBands,debug=True)
    plotInvvar('added-'+prefix,tractor)


    for i in range(itune):
        tractor.optimizeCatalogLoop(nsteps=1,srcs=[CG],sky=False)
        tractor.clearCache()
 #       saveBands('itune-%d-' % (i+1)+prefix,tractor,zr,flipBands,debug=True)
        plotInvvar('itune-%d-' % (i+1)+prefix,tractor)

    for i in range(ntune):
        tractor.optimizeCatalogLoop(nsteps=1,sky=True)

def test1():
	ps = PlotSequence('abell')

	run, camcol, field = 5115, 5, 151
	band = 'i'
	bands = [band]
	roi = (1048,2048, 0,1000)

	tim,tinf = st.get_tractor_image_dr9(run, camcol, field, band,
										roi=roi, nanomaggies=True)
	srcs = st.get_tractor_sources_dr9(run, camcol, field, band,
		roi=roi, nanomaggies=True,
		bands=bands)

	mags = [src.getBrightness().getMag(band) for src in srcs]
	I = np.argsort(mags)
	print('Brightest sources:')
	for i in I[:10]:
		print('  ', srcs[i])
	
	ima = dict(interpolation='nearest', origin='lower',
			   extent=roi)
	zr2 = tinf['sky'] + tinf['skysig'] * np.array([-3, 100])
	imb = ima.copy()
	imb.update(vmin=tim.zr[0], vmax=tim.zr[1])
	imc = ima.copy()
	#imc.update(vmin=zr2[0], vmax=zr2[1])

roi = [x0,x1,y0,y1]
    
    ra = 133.175
    dec = 33.4167
    itune = 10
    ntune = 5

    bands=['r','g','i','u','z']
    bandname = 'r'
    flipBands = ['r']

    rerun = 0

    TI = []
    TI.extend([st.get_tractor_image(run, camcol, field, bandname,roi=roi,useMags=True) for bandname in bands])
    sources = st.get_tractor_sources(run, camcol, field,bandname,roi=roi, bands=bands)

    timg,info = TI[0]
    photocal = timg.getPhotoCal()

    wcs = timg.getWcs()
    lvl = logging.DEBUG
    logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)

    tims = [timg for timg,tinf in TI]
    tractor = st.SDSSTractor(tims)
    tractor.addSources(sources)

    zr = np.array([-5.,+5.]) * info['skysig']

    print bands

band = 'r'
	roi = [1850, 1950, 620, 720]
	x0,y0 = roi[0],roi[2]

	# rd1 is near the position of the SDSS-catalog object
	# rd2 is four pixels above.
	rd1 = RaDecPos(225.67954, 11.265948)
	rd2 = RaDecPos(225.67991, 11.265852)

	lvl = logging.INFO
	logging.basicConfig(level=lvl, format='%(message)s', stream=sys.stdout)
	bandname = band

	sdssprefix = '%06i-%s%i-%04i' % (run, bandname, camcol, field)

	timg,info = st.get_tractor_image(run, camcol, field, bandname,
									 roi=roi)
	sources = st.get_tractor_sources(run, camcol, field, bandname,
									 roi=roi)

	wcs = timg.getWcs()
	for source in sources:
		x,y = wcs.positionToPixel(source, source.getPosition())
		print('  (%.2f, %.2f):' % (x+x0,y+y0), source)

	tractor = st.SDSSTractor([timg])
	tractor.addSources(sources)

	lnlorig = tractor.getLogLikelihood()
	zr = np.array([-5.,+20.]) * info['skysig']
	save(sdssprefix+'-orig', tractor, zr)

#   img.shape = data.shape
        #   del img.data
        #   del img.invvar
        #   R = super(Tractor2,self).getModelPatchNoCache(img, src)
        #   img.data, img.invvar = data,invvar
    
    

    #run,camcol,field = 7164,4,273
    #band='g'
    run,camcol,field = 2662, 4, 111
    band='i'
    roi=[0,300,0,300]
    im,info = st.get_tractor_image(run, camcol, field, band,
                                   useMags=True, roi=roi)
    sources = st.get_tractor_sources(run, camcol, field, band, roi=roi)
    tractor = Tractor2([im], sources)
    print tractor
    print tractor.getLogProb()
    tractor.freezeParam('images')

    p0 = tractor.getParams()
    tractor.setParams(p0)

    print
    print 'With Debug:'
    tractor.setParams(p0)
    tractor.mp = dmup
    t0 = Time()
    tractor.opt2()
    print 'With Debug:', Time()-t0
    print dpool.get_pickle_traffic_string()