#!/usr/bin/env python

from __future__ import absolute_import, division, print_function
from srcfinder_util import *

from skimage.filters.rank import median as medianfilt

import pylab as pl

NODATA=-9999

maxfetch = 1000 #np.inf
imeradinc = 5 # increment in meters to increase ime annuli radii
mergedists = [10,20,50]
mergediststr = str(mergedists)
use_src_loc = False

#srcxlsf='Source_Master_List_20170405_akt_tr_For_Bue.xlsx'
#sheetnames = ['Fall 2016_at_tr','High priority scenes']
#srcxlsf='Source_Master_List_Sort_20170515_bbedits.xlsx'
#srcxlsf='Source_Master_List_Sort_20170517-rmd_akt5.xlsx'
#sheetnames = ['Fall 2016_AVng']
#srcxlsf='Source_Master_List_Sort_20170614.xlsx'
#sheetnames = ['Fall 2016_AVng_Plumes']

colids = ['Line name','Candidate id',
          'Source id','Plume id',
          'Source Lat','Source Lon',
          'Plume Lat','Plume Lon']
colabv = ['lid','cid',
          'sid','pid',
          'slat','slon',
          'plat','plon']
colmap = dict(zip(colabv,colids))

configf = 'compute_ime.cfg'

def read_cfg(cfgf=configf):
    import ConfigParser

    config = ConfigParser.ConfigParser()
    config.read(cfgf)

    spreadsheet = config.get('Files', 'spreadsheet')
    worksheet = config.get('Files', 'worksheet')

    regen_all = config.get('Params', 'regen_all')
    regen_mtime = config.get('Params', 'regen_mtime')
    regen_lines = config.get('Params', 'regen_lines')

def write_cfg(config,cfgf=configf):
    if len(config.get('Params', 'regen_lines'))!=0:
        config.set('Params', 'regen_lines',[])
        
    with open(cfgf, 'w') as fid:
        config.write(fid)


def read_sources(srcxlsf,sheetname,filter=None,zone=None):
    from pandas import read_excel

    with open(srcxlsf) as fid:
        srcdata = read_excel(fid,sheetname=sheetname,skiprows=2)
    if srcxlsf == 'Source_Master_List_20170405_akt_tr_For_Bue.xlsx':
        with open(srcxlsf) as fid:
            srcdata = read_excel(fid,sheetname=sheetnames[1])
        srckeep = srcdata.ix[:,0]==1
    elif filter:
        filtercolid,filterval = filter.split(':')
        filtercol = np.array(srcdata.loc[:,filtercolid].values,dtype=str)
        filtercol = np.array(map(lambda s:s.strip(),filtercol),dtype=str)
                             
        srckeep=(filtercol==filterval.strip())
    else:
        srckeep = np.bool8(np.ones(len(srcdata)))
        
    # get column name mapping
    missing = []
    colloc = OrderedDict()
    print(srcdata.columns.values)
    for colid in colids:
        for colname in srcdata.columns.values:
            if colname.lower().startswith(colid.lower()):
                colloc[colid] = colname
                break
        if colid not in colloc:
            missing.append(colid)

    if len(missing)!=0:
        print('ERROR: missing columns in spreadsheet:',','.join(missing))
        sys.exit(0)
        
    srccols = colloc.values()        
    srcdata = srcdata.loc[srckeep,srccols]
    srcdata.columns = colloc.keys()    
    
    # print(srcdata)
    # utmcols = ['UTM X','UTM Y','UTM Zone','UTM Letter']
    # srclatlon = float32(srcdata[['Source Lat.','Source Lon.']])
    # srcutm = map(lambda latlon: latlon2utm(latlon[0],latlon[1],zone=zone),srclatlon)
    # plumelatlon = float32(srcdata[['Plume Lat.','Plume Lon.']])
    # plumeutm = map(lambda latlon: latlon2utm(latlon[0],latlon[1],zone=zone),plumelatlon)
    # for j,col in enumerate(utmcols):
    #     for i,v in enumerate(srcutm[j]):
    #         srcdata.ix[i,'Source '+col] = v
    #     for i,v in enumerate(plumeutm[j]):
    #         srcdata.ix[i,'Plume '+col] = v
    # print(srcdata)

    return srcdata

def ime(pixels_ppmm,ps):
    #           ppm(m)      ps^2        L/m^3       mole/L      kg/mole
    scalef = (1.0/1e6)*((ps*ps)/1.0)*(1000.0/1.0)*(1.0/22.4)*(0.01604/1.0)
    return (pixels_ppmm*scalef).sum()

def compute_ime(mfdat,mflab,lat,lon,mergedist,maxfetch,radps=imeradinc,
                doplot=False,outdir=None,nodata=NODATA):
    from scipy.spatial.distance import pdist,squareform

    assert(radps>0)
    mfmap = mfdat['mapinfo']
    ps = mfmap['xps']
    radinc = max(1,round(ps/radps))
    radincm = radps*radinc
    outdict = dict(detid=NODATA,detdim=NODATA,imevals=np.array([NODATA]),
                   detrgb=[],detmf=[],imergb=[],detpos=(NODATA,NODATA),                   
                   imemin=NODATA,imemax=NODATA,imesum=NODATA,
                   fetchm=NODATA,radincm=radincm,detbounds=[])

    outkeys = outdict.keys()

    j,i = np.int32(latlon2sl(lat,lon,mapinfo=mfmap))
    print('latlon2sl(',lat,lon,') ->',(j,i))
    plab = mflab[i,j]
    mfimg = mfdat['ch4mf']
    mfrgb = mfdat['rgb']
    mfmask = mfdat['nodata_mask']
    labmask = (mflab!=0) & (~mfmask)

    nodiam = 50 # extract a 100x100 context image by default for each position
    notpos = (i-nodiam,j-nodiam)
    notdim = (2*nodiam)+1
    notrgb = extract_tile(mfrgb,notpos,notdim,cval=0)*255
    detmf = extract_tile(mfimg,notpos,notdim,cval=0)
    outdict['detmf'] = detmf
    outdict['detrgb'] = np.uint8(notrgb)
    outdict['detpos'] = notpos
    outdict['detdim'] = notdim
    
    if plab==0 and labmask.any():
        #print('plab==0 for (i,j)=',(i,j))
        ij = np.c_[[i],[j]]
        ijdif = ij-np.c_[np.where(labmask)]
        ijoff = ijdif[np.argmin((ijdif**2).sum(axis=1))].squeeze()
        if ijoff.max()*ps > mergedist:
            return outdict

        i,j = i-ijoff[0],j-ijoff[1]
        plab = mflab[i,j]
        if plab!=plab or plab==0:
            return outdict
        
        #print('new (i,j)=',(i,j),'plab=',plab)

    pmask = mflab==plab    
    if not pmask.any():
        return outdict

    rc = np.where(pmask)
    # for tiny detections, fill pmask \pm 2
    nrc = len(rc[0])
    if nrc<=3:
        for idx in range(nrc):
            rci,rcj = rc[0][idx],rc[1][idx]
            r0=max(0,min(rci-2,mfimg.shape[0]-1))
            r1=max(0,min(rci+2,mfimg.shape[0]-1))+1
            c0=max(0,min(rcj-2,mfimg.shape[1]-1))
            c1=max(0,min(rcj+2,mfimg.shape[1]-1))+1            
            mflab[r0:r1,c0:c1]=plab*np.int32(mfimg[r0:r1,c0:c1]>0)
            
        pmask=mflab==plab
        rc = np.where(pmask)
        if len(rc[0])<=3:
            return outdict
        
    #print(radinc,j,i,rc)
    try:
        hrc = chull(np.c_[rc])
        #print('%d hull points'%len(hrc))
    except:
        print('qhull failed on %d input points'%len(rc[0]))
        raw_input()
        return outdict
    
    if 0:
        pl.ioff()
        fig,ax = pl.subplots(2,1,sharex=True,sharey=True,num=2)
        ax[0].imshow(mfimg.squeeze().transpose())
        ax[1].imshow(pmask.squeeze().transpose())
        ax[1].scatter(hrc[:,0],hrc[:,1],s=10,c='r',marker='o')
        pl.show()    

    pd = np.tril(squareform(pdist(hrc,metric='euclidean')))

    # compute fetch from distance equal to or slightly above pdmax (default=major axis)
    pdmax = min(round(maxfetch/ps),pd.max())
    pdmax = pd[(pd-pdmax)<=0].max()
    fetchm = pdmax*ps
    diam = int(np.ceil(pdmax))
    majaxi,majaxj = np.where(pd==pdmax)
    hmin,hmax = hrc[majaxi[0]],hrc[majaxj[0]]
    # print(diam)
    # print(hmin[0]-(i-diam),hmin[1]-(j-diam))
    # print(hmax[0]-(i-diam),hmax[1]-(j-diam))

    # extract pixel tiles
    tbuf = 5
    tdiam = int(max(diam+tbuf,nodiam))
    tpos = (i-tdiam,j-tdiam)
    tdim = (2*tdiam)+1
    tile = extract_tile(mfimg[...,np.newaxis],tpos,tdim,cval=0).squeeze()
    tmask = extract_tile(mfmask[...,np.newaxis],tpos,tdim,cval=0).squeeze()
    tlab = extract_tile(pmask[...,np.newaxis],tpos,tdim,cval=0).squeeze()    
    maxrad = bwdist(disk(2*tdiam))-tdiam
    maxrad = maxrad[tdiam:-tdiam,tdiam:-tdiam]
    tile[(tmask!=0) | (maxrad<=0)]=0
    circ = np.zeros(tile.shape)
    tlmask = (tlab!=0)

    radrange = range(1,tdiam+radinc+1,radinc)
    imevals = np.ones(len(radrange))*np.nan
    bar = progressbar('Computing IME for diam=%d'%diam,maxval=max(radrange))
    for i,rad in enumerate(bar(radrange)):
        imask = (maxrad>=rad) & (maxrad<rad+radinc) & (tlmask) 
        if not imask.any():
            imevals[i] = 0
            continue

        imeval = ime(tile[imask],ps)
        if imeval!=imeval: # ignore nans
            continue
        imevals[i] = imeval
        circ[imask] = imeval

    imevals = np.float64(imevals[~np.isnan(imevals)])

    cmap = 'YlOrRd'
    trgb = extract_tile(mfrgb,tpos,tdim,cval=0)*255
    tcmf = extract_tile(mfimg,tpos,tdim,cval=0)
    trgba = array2rgba(tile.reshape([tdim,tdim]),cmap=cmap,vmin=mfmin,vmax=mfmax)
    print('tlmask.sum(): "%s"'%str((tlmask.sum())))
    print('tile.shape: "%s"'%str((tile.shape)))
    print('trgb.shape: "%s"'%str((trgb.shape)))
    print('trgba.shape: "%s"'%str((trgba.shape)))
    ridx,cidx = np.where(tlmask)    
    trgb[ridx,cidx,:] = trgba[ridx,cidx,:-1]
    bridx,bcidx = np.where(outerboundaries(tlmask))
    detbounds = (bridx,bcidx)
    trgb[bridx,bcidx] = [0,255,0]
    trgb = np.uint8(trgb)            

    circ[circ==0] = np.nan
    crgb = array2rgba(circ.squeeze(),cmap=cmap,vmin=mfmin,vmax=mfmax)
    imemin,imemax,imesum = np.nanmin(circ),np.nanmax(circ),imevals.sum()
        
    outdict = dict(detid=plab,imevals=imevals,fetchm=fetchm,radincm=radincm,
                   detrgb=trgb,detmf=tcmf,imergb=crgb,detpos=tpos,detdim=tdim,
                   imemin=imemin,imemax=imemax,imesum=imesum,detbounds=detbounds)
    # sanity check
    for key in outkeys:
        if key not in outdict:
            print(key,'missing from return value')
            raw_input()
            
    return outdict


def mf_contours(mf,mfdetcc,mfmapinfo,mfmin=mfmin,levels=[500,1000,1500],
                mergedistm=50,doplot=False,return_sub=True):
    from skimage.segmentation import find_boundaries
    from skimage.measure import find_contours
    from geojson import MultiPolygon,Feature,FeatureCollection
    from matplotlib.colors import rgb2hex
    from scipy.ndimage.morphology import binary_fill_holes

    mfdist = int(round(mergedistm/mfmapinfo['xps']))
    # source_crs = dict(proj='longlat',ellps='WGS84',datum='WGS84',no_defs=True)
    # source_prop = dict(max_enhancement='float')
    # source_schema = dict(geometry='Polygon',properties=source_prop)
    
    # fill contours from min to max level (overwriting smaller levels)
    assert(mf.ndim == 2 or mf.shape[2]==1)
    assert(mfdetcc.ndim == 2 or mfdetcc.shape[2]==1)
    contourmap = np.zeros([mf.shape[0],mf.shape[1]],dtype=np.int32)
    np.set_printoptions(precision=9)
    mf_detids = np.unique(mfdetcc)[1:]
    print('mf_detids: "%s"'%str((mf_detids)))
    maxidx = np.zeros([3,len(mf_detids)],dtype=np.int32)
    levels_ordered = np.int32(sorted(levels))
    level_rgb = array2rgba(levels_ordered,cmap='YlOrRd')
    level_rgb = map(rgb2hex,level_rgb[:,:3]/255.0)
    mf_feats = []
    for detlab in mf_detids:
        detmask = mfdetcc==detlab
        (rmin,cmin),(rmax,cmax) = maskbbox(detmask,border=mfdist)
        mfsub =  mf[rmin:rmax,cmin:cmax].copy()
        masksub = detmask[rmin:rmax,cmin:cmax]
        detmaxv = mfsub[masksub].max()
        if detmaxv < mfmin:
            continue
        masksub = masksub*(mfsub>mfmin)
        if doplot:
            fig,ax = pl.subplots(1,2,sharex=True,sharey=True,num=1)
            ax[0].imshow(mfsub)
            ax[1].imshow(masksub)
            print(detlab,detmaxv,rmin,rmax,cmin,cmax)
            pl.show()
        
        ctrlonlat = []
        #print('cclab',cclab,'npts',ccmask.sum(),'max',detmaxv)
        detlevels = levels_ordered[levels_ordered<=detmaxv]
        selem = disk(3)
        nlevel = 3
        for li,leveli in enumerate(detlevels):        
            detlevelmask = binary_fill_holes(masksub & (mfsub>=leveli))
            detlevelmask = medianfilt(bwdilate(detlevelmask,selem=selem),selem)
            detlevelmask = np.float32(detlevelmask!=0)
            #pl.imshow(detlevelmask)
            #pl.show()
            detlevelr,detlevelc = np.where(detlevelmask!=0)
            contourmap[detlevelr+rmin,detlevelc+cmin] = leveli
            for ilevel in range(1,nlevel+1):
                flevel = ilevel/float(nlevel+1)
                cclevelctr = find_contours(detlevelmask,flevel,
                                           fully_connected='low',
                                           positive_orientation='low')
                nctr = len(cclevelctr)
                if nctr==0:
                    print('no contours found, trying alternate connectivity')
                    cclevelctr = find_contours(detlevelmask,flevel,
                                               fully_connected='high',
                                               positive_orientation='high')
                    nctr = len(cclevelctr)
                    
                if nctr!=0:
                    # found contours
                    print('%d contours found at level %.3f'%(nctr,flevel))
                    break
            print('contours for id=',detlab,'level=',leveli,'# contours=',len(cclevelctr),
                  'total area=',np.count_nonzero(detlevelmask))
            for ctr in cclevelctr:
                ctrlat,ctrlon = sl2latlon(cmin+np.round(ctr[:,1]),
                                          rmin+np.round(ctr[:,0]),
                                          mapinfo=mfmapinfo)
                ctrlonlat.append(zip(ctrlon,ctrlat))
        ccmaxrow,ccmaxcol = np.where(masksub & (mfsub==detmaxv))
        contourmap[rmin+ccmaxrow,cmin+ccmaxcol] = int(detmaxv)
        ccmaxlat,ccmaxlon = sl2latlon(rmin+ccmaxrow,cmin+ccmaxcol,#rmin+ccmaxcol,cmin+ccmaxrow,
                                      mapinfo=mfmapinfo)
        ccprop = dict(levels=['%d'%l for l in detlevels],
                      maxenhance=str((detmaxv,[ccmaxlat[0],ccmaxlon[0]])))
        ccfeat = Feature(properties=ccprop,geometry=MultiPolygon(ctrlonlat))
        mf_feats.append(ccfeat)
        if doplot:
            fig,ax = pl.subplots(1,2,sharex=True,sharey=True,num=1)
            ax[0].imshow(mfdetcc[rmin:rmax,cmin:cmax])
            ax[1].imshow(contourmap[rmin:rmax,cmin:cmax])
            pl.suptitle('detection id=%d'%detlab)
            pl.show()
    for i,poly in enumerate(mf_feats):
        print('mf_feats polygon',i,len(poly),'\n')
    mf_feats = FeatureCollection(mf_feats)    
    mf_feats = FeatureCollection(mf_feats)
    if return_sub:
        (rmin,cmin),(rmax,cmax) = maskbbox((contourmap!=0),border=mfdist)
        return contourmap[rmin:rmax,cmin:cmax],mf_feats
    return contourmap, mf_feats

def dtstr():
    import datetime as dtime
    _now = dtime.datetime.now()
    return _now.strftime('%a %B %d, %H:%M:%S %Y')

def process_mfimg(mfinfile,outdir,srcdata,colmap,args,
                  ctr_smooth=0,NODATA=NODATA):
    savefilt = True #args.savefilt
    mergedists = args.mergedists
    overwrite = args.clobber
    doplot = True #args.plot
    maxfetch = args.fetchmax
    nodetfilt = args.nodetfilt

    if not pathexists(outdir):
        mkdir(outdir)

    lineids = srcdata[colmap['lid']]

    imgid = filename2flightid(mfinfile)
    imgdate = filename2flightdate(mfinfile)    
    srcplumes = srcdata[lineids==imgid]

    mfdir,mffile = pathsplit(mfinfile)

    statsf = pathjoin(outdir,mffile+'_ime.txt')
    if pathexists(statsf) and not overwrite:
        print(statsf,'already exists')
        return 
        
    #detfiltf = args.detfilt
    #savefilt = False if (detfiltf and pathexists(detfiltf)) else savefilt

    nsrc = len(srcplumes)
    print('%d plumes for image %s'%(nsrc,imgid))
    if nsrc == 0:
        return

    mfdata = loadmfdata(mfinfile)
    ch4mf = mfdata['ch4mf']
    ch4rgb = mfdata['rgb']
    mfmask = mfdata['nodata_mask']
    mfmap = mfdata['mapinfo']
    mf_ps = mfmap['xps']
    mfignore = np.where(mfmask)

    # if detfiltf:
    #     detdata = loaddetfilt(detfiltf.strip())
    #     ch4det = detdata['ch4det']
    #     detcc = imlabel((ch4det!=0))
    # else:
    det_outf = detkde_outf = detcomp_outf = None
    if savefilt:
        out_suf = '_filt_det_%d_%d'%(mfmin,mfmax)
        outkde_suf = '_filt_det_kde_%d_%d'%(mfmin,mfmax)
        outcomp_suf = '_filt_det_ccomp_%d_%d'%(mfmin,mfmax)
        out_base = splitext(pathsplit(mfinfile)[1])[0]
        det_outf = pathjoin(outdir,out_base+out_suf)
        detkde_outf = pathjoin(outdir,out_base+outkde_suf)
        detcomp_outf = pathjoin(outdir,out_base+outcomp_suf)

    skip_kde = False
    if nodetfilt:
        areamin = 0
    else:
        areamin = minarea

    ch4det,detcc = filtdet(ch4mf,mfmask,mfmap,minarea=areamin,
                           skip_kde=skip_kde,
                           det_outf=det_outf,
                           kde_outf=detkde_outf,
                           ccomp_outf=detcomp_outf)
        
    #detids = ccomp2detids(detcc,mf_ps)
    #point_ids = detids['point_ids']
    #diffuse_ids = detids['diffuse_ids']

    mfmapstr = mapdict2str(mfmap)
    if 0:
        from skimage.future import graph
        ch4cc = sliczero(ch4mf[:,:,np.newaxis],areamin**2)
        g = graph.rag_mean_color(ch4mf, ch4cc, mode='similarity')
        gids = graph.cut_normalized(ch4cc, g)
        print(gids.min(),gids.max())
        if savefilt:
            gids_outf = det_outf+'_idsncut'
            array2img(gids_outf,gids,mapinfostr=mfmapstr,overwrite=True)
            print('saved',gids_outf), raw_input()
        

    det_ids = {}
    for i,dist in enumerate(sorted(mergedists)):
        print('computing detection map for mergedist=%6.3f'%dist)        
        if i>0:
            # use previously merged ids computed using smaller threshold 
            dprev=mergedists[i-1]
            ids = np.int32(mergelabels(det_ids[dprev],int(round(dist/mf_ps))))
            if (det_ids[dprev]==ids).all():
                warn('detection maps for mergedist=%7.3f and %7.3f equal'%(dprev,dist))
        else:
            ids = np.int32(mergelabels(detcc,int(round(dist/mf_ps))))

        srcdir = pathjoin(outdir,'mergedist%d'%int(dist))
        if not pathexists(srcdir):
            print('creating directory',srcdir)
            mkdir(srcdir)
            
        if savefilt:
            det_idf = pathsplit(det_outf)[1]+'_ids%d'%dist
            ids_outf = pathjoin(srcdir,det_idf)
            array2img(ids_outf,ids,mapinfostr=mfmapstr,overwrite=True)
        det_ids[dist] = ids
            
    sids = srcplumes[colmap['sid']].values
    pids = srcplumes[colmap['pid']].values
    cids = srcplumes[colmap['cid']].values
    slats = srcplumes[colmap['slat']].values
    slons = srcplumes[colmap['slon']].values
    plats = srcplumes[colmap['plat']].values
    plons = srcplumes[colmap['plon']].values

    imehdrs = ['IME%d (kg)'%d for d in mergedists]
    fetchhdrs = ['Fetch%d (m)'%d for d in mergedists]
    detidhdrs = ['DetId%d'%d for d in mergedists]
    hdrvals = [colmap['cid'],colmap['sid'],colmap['pid']]+imehdrs+fetchhdrs+detidhdrs
    lines = ['# '+', '.join(hdrvals)]
    mf_ctr = np.zeros_like(ch4mf)
    for i in range(nsrc):
        sid,cid,pid = sids[i],cids[i],pids[i]
        
        plat,plon = plats[i],plons[i]
        slat,slon = slats[i],slons[i]

        lat,lon = (slat,slon) if use_src_loc else (plat,plon)           
        
        print('source',i,sid,cid,pid,
              '\nsource (lat,lon)=',(slat,slon),
              '\nplume (lat,lon)=',(plat,plon),
              '\nusing source location=',use_src_loc)
        imesumi,fetchmi,detidi = [],[],[]
        for dist in mergedists:
            idmask = det_ids[dist]
            print('mergedist=',dist,'detection count=',idmask.max())
            imeout = compute_ime(mfdata,idmask,lat,lon,dist,maxfetch)
            
            srcdir = pathjoin(outdir,'mergedist%d'%int(dist))
            detpos, detdim = imeout['detpos'],imeout['detdim']
            srcpos = 'r%d_c%d'%detpos
            srcbase = pathjoin(srcdir,'_'.join([imgid,sid,srcpos]))
            rgbimgf = srcbase+'_rgb.png'
            qlimgf = srcbase+'_rgbql%d_%d.png'%(mfmin,mfmax)
            ppmimgf = srcbase+'_cmf%d_%d.png'%(mfmin,mfmax)
            imeimgf = srcbase+'_ime.png'
            ctrimgf = srcbase+'_ctr.png'
            ctrgtif = srcbase+'_ctr.tif'
            ctroutlimgf = srcbase+'_ctr_outl.png'            
            figimgf = srcbase+'_fig.pdf'
            ctrf = srcbase+'_ctr.geojson'

            detrgb = imeout['detrgb']
            imsave(rgbimgf,detrgb,verbose=True)

            detmf = imeout['detmf']
            detmask = (detmf>=mfmin)
            print(detmf.shape,detrgb.shape)
            detql = rgbdet2ql(detrgb,detmf,detmask)
            imsave(qlimgf,detql,verbose=True)

            ppmql = detql.copy()
            detidx = np.where(~detmask)
            if len(detidx[0])>0:
                igr,igc = detidx[0],detidx[1]
                ppmql[igr,igc,:] = 0
            imsave(ppmimgf,ppmql,verbose=True)
            
            detid = imeout['detid']
            if detid==NODATA:
                msg = 'WARNING (sourceid=%s, plumeid=%s):'%(sid,pid)
                msg+= 'no pixels >=%7.2fppmm within %dm of (%f,%f), '%(mfmin,dist,
                                                                       lat,lon)

                msg+='observed ppmm range: %s, '%(str(extrema(detmf)))
                msg+='observed [%d,%d] area: %d'%(mfmin,mfmax,(detmf>mfmin).sum())
                print(msg)                
                continue

            
            imergb = imeout['imergb']
            imsave(imeimgf,imergb,verbose=True)
            if doplot:
                fig,ax = pl.subplots(1,2,sharex=True,sharey=True,num=1)
                ax[0].imshow(imeout['detrgb'])
                ax[0].set_xlabel('CH$_{4}$ MF ppmm Min=%-6.3f, Max=%-6.3f'%(mfmin,mfmax))
                ax[1].imshow(imeout['imergb'])
                ax[1].set_xlabel('IME $\Delta$=%4.2fm Min=%-6.3f, Max=%-6.3f'%(imeout['radincm'],
                                                                               imeout['imemin'],
                                                                               imeout['imemax']))
                pl.suptitle('IME=%-6.3f kg, Fetch=%7.2fm'%(imeout['imesum'],imeout['fetchm']))
                pl.setp([a.get_xticklabels() for a in fig.axes], visible=False)
                pl.savefig(figimgf)
                pl.clf()
                print('saved',figimgf)

            detmask = np.int32(idmask==imeout['detid'])
            mf_ctr, mf_feats = mf_contours(ch4mf,detmask,mfmap,mfmin=mfmin,
                                           levels=[500,1000,1500],
                                           return_sub=False)
            
            with open(ctrf,'w') as fid:
                print(str(mf_feats),file=fid)
            ctrvalid = extract_tile((~mfmask) & (mf_ctr==0),detpos,detdim,cval=0).squeeze()
            ctrvalid = np.where(ctrvalid)
            ctrmask = extract_tile(mfmask,detpos,detdim,cval=0).squeeze()
            #mf_ctr = mf_ctr.squeeze()
            #ctr_shape = list(mf_ctr.shape)
            #print('mf_ctr.shape: "%s"'%str((mf_ctr.shape)))
            if 0:
                mfignore = np.where(mfmask)
                mfvalid = np.where(ctrvalid)
                ctr_rgba = array2rgba(mf_ctr,cmap='YlOrRd',vmin=500,vmax=1500) #.reshape(ctr_shape+[4])
                #for i in range(ctr_smooth):
                #    for j in range(3):
                #        ctr_rgba[...,j] = medianfilt(ctr_rgba[...,j], disk(3))
                
                ctr_rgba[mfignore[0],mfignore[1],-1] = 0
                ctr_rgba[mfvalid[0],mfvalid[1],:-1] = 0 # black bg
                ctr_rgba[mfvalid[0],mfvalid[1],-1] = 0 # transparent bg            
                ctr_rgba = extract_tile(ctr_rgba,detpos,detdim,cval=0)
                
            detctr = extract_tile(mf_ctr,detpos,detdim,cval=0).squeeze()
            detctr[ctrmask]=0
            ctr_rgba = array2rgba(detctr,cmap='YlOrRd',vmin=500,vmax=1500) #.reshape(ctr_shape+[4])
            ctrignore = np.where(ctrmask)
            # transparent bg
            ctr_rgba[ctrignore[0],ctrignore[1],-1] = 0
            ctr_rgba[ctrvalid[0],ctrvalid[1],-1] = 0 
            print('ctr_rgba.shape: "%s"'%str((ctr_rgba.shape)))
            imsave(ctrimgf,ctr_rgba,verbose=True)
            tile2geotiff(ctr_rgba,detpos,ctrgtif,mfinfile)
            #outl_rgba = ctr_rgba.copy()
            # rgbctr = extract_tile(detrgb,detpos,detdim,cval=0)
            # outl_rgb = rgbdet2ql(rgbctr,detctr,ctrmask)
            
            # oidx = np.where(outerboundaries(np.uint8(detctr.any(axis=2))))
            # if len(oidx)!=0:
            #     outl_rgb[oidx[0],oidx[1]] = [0,255,0]
            #     imsave(ctroutlimgf,outl_rgb,verbose=True)
            imesumi.append(imeout['imesum'])
            fetchmi.append(imeout['fetchm'])
            detidi.append(imeout['detid'])
            
        outvals = map(str,[cid,sid,pid]+imesumi+fetchmi+detidi)
        lines.append(', '.join(outvals))

        
    with open(statsf,'w') as fid:
        print('\n'.join(lines),file=fid)
    print('wrote',statsf)

if __name__ == '__main__':
    import argparse

    parser = argparse.ArgumentParser(description="CMF Prediction Summary")
    
    #parser.add_argument('-d','--detfilt', type=str,  
    #                    help='path to filtered detection image (computed from mfinput by default)')

    mergedisthelp='List of one or more distances (in meters) to merge neighboring blobs (default=%s)'%mergediststr
    parser.add_argument('-m','--mergedists', nargs='+', type=float,
                        help=mergedisthelp, default=mergedists)
    
    #parser.add_argument('-p','--plot', action='store_true',
    #                    help='Save IME plots to outdir')

    parser.add_argument('-c','--clobber', action='store_true',
                        help='Clobber existing outputs')      

    #parser.add_argument('-s','--savefilt', action='store_true', 
    #                    help='Save filtered detection image to outdir')      

    parser.add_argument('-o','--outdir', type=str, default='./ime/',
                        help='Save filtered detection image to OUTDIR')

    parser.add_argument('-f','--fetchmax', type=float, default=maxfetch,
                        help='Maximum fetch value (default=unlimited)')

    parser.add_argument('--nodetfilt', action='store_true',
                        help='Do not filter out small or weakly connectedd detections')
    
    parser.add_argument('-i','--inputfile', type=str,  metavar='inputfile',
                        help='CMF image')

    parser.add_argument('-l','--listfile', type=str,  metavar='listfile',
                        help='List of CMF images to process')    

    parser.add_argument('--randomize', action='store_true', 
                        help='Randomize sequence processing order')      
    
    parser.add_argument('--filter', type=str,  metavar='filter',
                        help='only compute IME for spreadsheet rows matching filter')    
    
    parser.add_argument('xlsfile', type=str,  metavar='xlsfile',
                        help='Source master list spreadsheet')

    parser.add_argument('sheetname', type=str,  metavar='sheetname',
                        help='Source master list spreadsheet sheet name')    

    args = parser.parse_args()
    outdir = args.outdir
    inputfile = args.inputfile
    overwrite = args.clobber
    listfile = args.listfile
    xlsfile = args.xlsfile
    sheetname = args.sheetname
    srcfilter = args.filter
    randomize = args.randomize
    
    srcdata = read_sources(xlsfile,sheetname,filter=srcfilter)
    lineids = srcdata[colmap['lid']]
    
    # remove nans
    keepmask = lineids==lineids 
    srcdata = srcdata[keepmask]
    lineids = np.unique(lineids[keepmask])

    if randomize:
        rperm = np.random.permutation(len(lineids))
        lineids = lineids[rperm]
    
    if inputfile:
        mfinfiles = [inputfile]
    elif listfile:
        mfinfiles = np.loadtxt(listfile,dtype=str)
    else:
        from socket import gethostname as hostname
        if hostname().startswith('valis'):
            cmfdir = './data'
        else:
            cmfdir = '/lustre/ang'
        cmfdir = pathjoin(cmfdir,'y16/cmf/ort')
        cmfsuf = '_cmf_v1n2_img'
        mfinfiles = [pathjoin(cmfdir,lid+cmfsuf)
                     for lid in lineids]        

    nfiles = len(mfinfiles)
    for i,mfinfile in enumerate(mfinfiles):        
        if 1: #try:
            if not pathexists(mfinfile):
                print(mfinfile,'not found')
                continue
            imgid = filename2flightid(mfinfile)
            mfoutdir = pathjoin(outdir,imgid)
            if randomize and pathexists(mfoutdir) and not overwrite:
                # only skip directory paths when randomize is on,
                # otherwise regenerate outputs whenever _ime.txt is missing
                print('path',mfoutdir,'exists, skipping')
                continue
            print(dtstr()+': started processing mfinfile=',mfinfile,'(%d of %d)'%(i+1,nfiles))
            process_mfimg(mfinfile,mfoutdir,srcdata,colmap,args)
            print('-'*65)
            print('')
        #except Exception as e:
        #    print('EXCEPTION:',str(e))