Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mrao.cam.ac.uk/~bn204/alma/memo-infer-587/sma_retrieval.py Дата изменения: Fri Aug 10 20:17:35 2012 Дата индексирования: Tue Oct 2 12:39:24 2012 Кодировка: Поисковые слова: чоеыойе рмбоефщ

"""
Bojan Nikolic ,
Initial version December 2008
Revised 2009

Retrieval of atmospheric parameters from SMA data
"""

import math
import os
import numpy
import pyfits

import pickle

import iofits4

import smadatasets
import sma_plots
import air
from sma_skydip import _dip1_times
import bnmin1io
from likelihood_bplot import PlotMCMC

def stdRadiometer(hin,
dickeAB="A",
mask=None):
"""
Convert radiometer data to standard format
"""
if mask is None:
mask=numpy.ones(len(hin.data),
numpy.bool)
if dickeAB=="A":
d=numpy.array([hin.data.field("ch%i"%i)[mask]for i in range(4)],
numpy.float64)
elif dickeAB=="B":
d=numpy.array([hin.data.field("ch%i"%i)[mask]for i in range(4,8)],
numpy.float64)
d=d[::-1,:]
else:
raise "unknown radiometer "+dickeAB
return d

def getRadioData(dset,
dickeAB="A"):
"""
Get radiometer data in suitable format

:param dickeAB: Select the "A" or the "B" radiometers within the
prototype unit

:returns: Tuple (time, elevation, data)

"""
fname=dset[1]
hin=pyfits.open(fname)[1]
el=hin.data.field("el")
time=hin.data.field("time")

mask=pyfits.open(fname[:-5]+".mask.fits")[1].data.field("SCal")==1

el=el[mask]
time=time[mask]
d=stdRadiometer(hin,
mask=mask,
dickeAB=dickeAB)
return time, el, d

def getDipData(fnamein,
times,
dickeAB="A",
delta=0.0002):
"""
Extract dip data, i.e., clusters of points around the times
parameter

example:
>>> t,e,d=getDipData("o/temp/skydip1.fits", _dip1_times)
>>> x2=dipRetrieve(e, d,nsample=800000,fitTTerm=True)

"""
hin=pyfits.open(fnamein)[1]
din=hin.data
tin=din.field("time")
elin= din.field("el")
mask=numpy.zeros(len(tin),
numpy.bool)
for t in times:
mask=numpy.logical_or(mask,
numpy.logical_and(tin>t-delta,
tin d=stdRadiometer(hin,
mask=mask,
dickeAB=dickeAB)
d=d[::-1]
return (tin[mask],
elin[mask],
d.transpose())



def retAtmPars(dset,
dtp=100,
ns_base=10000,
dickeAB="A",
trange=None,
el_override=None):
"""
Retrieve atmospheric parameters from SMA campaign observations.

:param dset: The data set to analyse

:param dtp: Number of samples to skip between each retrieval (set
to 1 to get a retrieval for each sample). If dtp is is a list,
then retrieve for these time points only

:param ns_base: Base scaling of the number of proposals in the
MCMC to make

:param el_override: Ignore the elevation record in data and use
the value of this parameter instead

Example:

>>> retAtmPars(smadatasets.sfeb17)

>>> x=retAtmPars(smadatasets.sfeb17, dtp=1000, trange=(17.5,17.55), ns_base=20000)

"""
t,el,d=getRadioData(dset,
dickeAB=dickeAB)
if trange:
mask=numpy.logical_and(t>trange[0],
t el=el[mask]
d=d[:,mask]
n=len(el)

res=[]

if type(dtp) is list:
tpl=dtp
else:
tpl=range(0,n,dtp)

for tp in tpl:
dp=d[::-1,tp]
cel=el[tp]
if el_override is not None:
cel=el_override
print cel,dp
if len(res) is 0:
nsample=ns_base*10
ic=[1.0, 270, 600]
else:
nsample=ns_base
ic=res[-1][1][-1]
r=air.AbsRetrieve(d=dp,
ic=ic,
nsample=nsample,
prior=air.MK_flat_pr,
coupling=0.91,
TSpill=265,
el=cel)
res.append((cel, r))
return res

def retAtmSave(dset,
fnameout,
**args):
r=retAtmPars(dset,
**args)
for i, (el, chain) in enumerate(r):
bnmin1io.fSaveChain(["n", "T", "P"],
chain,
fnameout+("-%04i.fits" % i),
kwrds={"el" : el,
"proc000" : "Atmospheric parameter retrieval"})


def retSimul(ic=[1.0, 270, 580],
nsample=100000,
prior=None):
"""
Retrieve from a simulated absolute temperature observation and
save to FITS file
"""
fnameout="cache/simul"
if prior is not None:
if prior is "tight":
prior=air.MK_tight_pr
fnameout+="MKTight"
else:
prior=air.MK_flat_pr
fnameout+="MKPrior"
else:
prior=air.MK_weak_pr
r=air.AbsRetrieve(d=None,
ic=ic,
nsample=nsample,
prior=prior,
PSigma=0.5,
TSigma=0.2,
)
bnmin1io.fSaveChain(["n", "T", "P"],
r,
fnameout+".fits")

def plotSimRet(prior=False):
fnamein="cache/simul"
pref="abs"
if prior:
if prior is "tight":
fnamein+="MKTight"
pref+="MKTight"
else:
fnamein+="MKPrior"
pref+="MKPrior"
pl=["n", "T", "P"]
d=loadChain(fnamein+".fits",
pl)
PlotMCMC(d,pl,
pref=pref,
d2bins=100)

pl2=[("dLdT%i"%i) for i in range(1,5)]
d=loadChain(fnamein+"-dtdl.fits",
pl2)
PlotMCMC(d,
pl2,
pref=pref,
d2bins=100)

def plotSingle(infix=""):
fnamein="cache/f17Single%s-0000" % infix
pref=("f17mid"+infix)
pl=["n", "T", "P"]
d=loadChain(fnamein+".fits",
pl)
PlotMCMC(d,pl,
pref=pref,
d2bins=100,
burn=30000)

pl2=[("dLdT%i"%i) for i in range(1,5)]
d=loadChain(fnamein+"-dtdl.fits",
pl2)
PlotMCMC(d,
pl2,
pref=pref,
d2bins=100)





def dTdLSave(fnamein,
fnameout,
coupling=0.91,
params=["n", "T", "P"],
trunc=None):
"""
Turn a FITS file chain of model parameters to a FITS file of dTdLs
"""
fin=pyfits.open(fnamein)
if fin[0].header.has_key("el"):
el=fin[0].header["el"]
else:
el=math.pi*0.5
d=loadChain(fnamein,
params)
dtdl=air.calcdTdL(d,
coupling=coupling,
el=el)
coldefs = [
pyfits.Column("dLdT%i"%i,
"E",
"mm/K",
array=1.0/dtdl[:,i-1]) for i in range(1,5)]
tabout=pyfits.new_table(coldefs)
fout=iofits4.PrepFitsOut("")
fout.append(tabout)
iofits4.Write(fout,
fnameout,
overwrite=1)
return dtdl



def removeBurn(r,
f=0.5):
"""
Remove the burn-in part of a set of retrievals
"""
rr=[]
for x in r:
n=len(x[1])
rr.append((x[0],x[1][f*n:]))
return rr

def globalMax(r,i):
return max([x[1][:,i].max() for x in r])

def globalMin(r,i):
return min([x[1][:,i].min() for x in r])


def globalHist(r,
i,
nbins=100):
"""
Compute the sequence of histogram for full observation set, on
same bins

Example:
>>> m=globalHist(rr,2,nbins=200)
>>> pylab.matshow(m.transpose())
"""
res=[]
range=(globalMin(r,i),globalMax(r,i))
for x in r:
h,bins=numpy.histogram(x[1][:,i],
range=range,
bins=nbins)
res.append(h/float(h.sum()))
return numpy.array(res), bins

def loadChain(fnamein,
params):
fin=pyfits.open(fnamein)
d=numpy.array([fin[1].data.field(x) for x in params],
numpy.float64)
d=d.transpose()
return d

def loadRetSeq(pref,
params):
"""
Load a sequence of chains saved as FITS files

Example:

>>> r=loadRetSeq("cache/feb17-za-A.p-", ["n", "T", "P"] )

"""
i=0
res=[]
while True:
fnamein=pref+ ("%04i.fits"%i)
if os.access(fnamein, os.F_OK):
fin=pyfits.open(fnamein)
el=float(fin[0].header["el"])
d=loadChain(fnamein, params)
res.append([el,d])
i=i+1
else:
break
return res


def plotRetSeq(r,
params,
pref="",
):
"""
Plot results of sequential retrieval parameters

Example:

>>> plotRetSeq(pickle.load(open("cache/feb17-3modret-B.p")) , ["n", "T", "P"],pref="simpleB")

>>> plotRetSeq(pickle.load(open("cache/test-feb17-3modret-dtdl-A.p")) , ["dTdL1", "dTdL2", "dTdL3","dTdL4"],pref="test")

"""
rr=removeBurn(r)
for i,p in enumerate(params):
m,bins=globalHist(rr,i,nbins=200)
sma_plots.plotParRetrSeq(m,
bins,
p,
"o/%s-%s-seq.eps" % (pref, p ))



def doRetrieval(dset,
pref,
trange=None,
doParam=True,
doDTDL=True,
dtp=25,
ns_base=10000):
"""
Helper routine which carries out retrieval on an SMA data set and
writes out the results.

:param doParam: Retrieve the model parameters?

"""
for dickeAB in ["A", "B"]:
if doParam:
retAtmSave(dset,
pref+"-%s.p" % dickeAB,
trange=trange,
dtp=dtp,
dickeAB=dickeAB,
ns_base=ns_base)
if doDTDL:
dTdLSeqSave(pref+"-%s.p" % dickeAB,
pref+"-dtdl-%s.p" % dickeAB)

def analyseABComparison():
"""
Do a few points with very long chains to check difference between
A and B receivers of the Dicke radiometer
"""
dtp=numpy.arange(250,270,5)*25
doRetrieval(smadatasets.sfeb17,
"cache/ABComparison",
dtp=list(dtp),
ns_base=100000)


def Analyse(doall=False):

if False or doall:
retSimul(nsample=2000000)
retSimul(nsample=2000000,
prior=True)
x=dTdLSave("cache/simul.fits",
"cache/simul-dtdl.fits",
coupling=1.0)
x=dTdLSave("cache/simulMKPrior.fits",
"cache/simulMKPrior-dtdl.fits",
coupling=1.0)

if False or doall:
retSimul(nsample=2000000, prior="tight")
dTdLSave("cache/simulMKTight.fits",
"cache/simulMKTight-dtdl.fits",
coupling=1.0)

if False or doall:
plotSimRet()
plotSimRet(prior=True)

if True or doall:
retAtmSave(smadatasets.sfeb17,
"cache/f17Single",
dtp=2000,
trange=(17.5,17.55),
ns_base=200000)
dTdLSave("cache/f17Single-0000.fits",
"cache/f17Single-0000-dtdl.fits",
coupling=0.91)

retAtmSave(smadatasets.sfeb17,
"cache/f17Single-beg",
dtp=2000,
trange=(17.1,17.15),
ns_base=200000)
dTdLSave("cache/f17Single-beg-0000.fits",
"cache/f17Single-beg-0000-dtdl.fits",
coupling=0.91)

retAtmSave(smadatasets.sfeb17,
"cache/f17Single-end",
dtp=2000,
trange=(17.8,17.85),
ns_base=200000)
dTdLSave("cache/f17Single-end-0000.fits",
"cache/f17Single-end-0000-dtdl.fits",
coupling=0.91)

plotSingle()
if False or doall:
doRetrieval(smadatasets.sfeb17,
"cache/feb17-za",
trange=(17.1,17.85),
ns_base=50000,
doDTDL=False,
dtp=10)