#!/usr/bin/env python from __future__ import absolute_import, division, print_function import matplotlib from socket import gethostname # use agg backend unless we're running on one of the hosts below hostname = gethostname() if hostname not in ('valis','gibson'): matplotlib.use('agg') from srcfinder_util import * from skimage.filters.rank import median as medianfilt import pylab as pl from glob import glob NODATA=-9999 MINRAD=75 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'] # define paths and filename regexps for each platform (ang,avcl,etc.) if hostname.startswith('valis'): cmfdir_ang = pathjoin(expanduser('~/Research/srcfinder'),'data') cmfdir_avcl = pathjoin(expanduser('~/Research/srcfinder'),'data_avcl') else: cmfdir_ang = '/lustre/ang' cmfdir_avcl = '/lustre/athorpe/ch4/compute_ime_new_test/AVIRIS_Scenes' calids_ang = ['v1f','v1n2'] cmfsuf_ang = '_cmf_*_img' # replace '*' with calids_ang calids_avcl = ['b_sc01'] cmfsuf_avcl = '_*_ort_img_cmfv' # replace '*' with calids_ang # search calids for ang+avcl by default calids = calids_ang+calids_avcl # override srcfinder_util minarea minarea_ang = 9 minarea_avcl = 3 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)) colorf = 'levels_rgba.txt' def parse_colors(colorf): colors = np.loadtxt(colorf,delimiter=':',dtype=str,comments='#') levels_ppmm = np.float32(colors[:,0]) levels_rgba = map(lambda s: s.split(','),colors[:,1]) levels_rgba = np.float32(levels_rgba)/255 return levels_ppmm,levels_rgba # 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,filterstr=None,zone=None): from pandas import read_excel from collections import OrderedDict with open(srcxlsf) as fid: srcdata = read_excel(fid,sheet_name=sheetname,skiprows=2) if filterstr: filtercolid,filterval = filterstr.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() 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_scale(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 scalef def ime(pixels_ppmm,ps): return pixels_ppmm.sum()*ime_scale(ps) def compute_ime(mfdat,mflab,lat,lon,mergedist,maxfetch,mfmin=mfmin, 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 = ps*radinc outdict = dict(detid=NODATA,detdim=NODATA,detpos=(NODATA,NODATA), imevals=np.array([NODATA]),cmfnums=np.array([NODATA]), cmfvals=np.array([NODATA]),imerads=np.array([NODATA]), detrgb=[],detql=[],detmf=[],imergb=[], 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['image'] mfrgb = mfdat['rgb'] mfmask = mfdat['nodata_mask'] labmask = (mflab!=0) & (~mfmask) & (mfimg>=mfmin) nodiam = MINRAD # extract a 2MINRAD+1 x 2MINRAD+1 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) detrgb = np.uint8(notrgb) outdict['detmf'] = detmf outdict['detrgb'] = detrgb outdict['detql'] = detrgb 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: bufrad=2 for idx in range(nrc): rci,rcj = rc[0][idx],rc[1][idx] r0=max(0,min(rci-bufrad,mfimg.shape[0]-1)) r1=max(0,min(rci+bufrad,mfimg.shape[0]-1))+1 c0=max(0,min(rcj-bufrad,mfimg.shape[1]-1)) c1=max(0,min(rcj+bufrad,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) nrc = len(rc[0]) if nrc<=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'%nrc) if nrc < 10: print(rc) 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 = np.arange(1,tdiam+radinc+1,radinc) imevals = np.ones(len(radrange))*np.nan cmfvals = np.zeros(len(radrange)) cmfnums = np.zeros(len(radrange)) bar = progressbar('Computing IME for diam=%d'%diam,maxval=max(radrange)) for i,rad in enumerate(bar(radrange)): imask = (maxrad>=rad) & (maxrad0) imevals = np.float64(imevals[imekeep]) imerads = radrange[imekeep]/ps cmfnums = cmfnums[imekeep] cmfvals = cmfvals[imekeep] 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('total pixels used in IME calculation: "%s"'%str((tlmask.sum()))) print('output tile shape: "%s"'%str((trgba.shape))) ridx,cidx = np.where(tlmask) trgb[ridx,cidx,:] = trgba[ridx,cidx,:-1] # draw a green border surrounding the plume extent? bridx,bcidx = np.where(outerboundaries(bwdilate(tlmask,selem=disk(2)))) detbounds = (bridx,bcidx) draw_bounds = False if draw_bounds: trgb[bridx,bcidx] = [0,255,0] detql = 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,cmfnums=cmfnums,imerads=imerads, cmfvals=cmfvals,fetchm=fetchm,radincm=radincm,detpos=tpos, detdim=tdim,detrgb=detrgb,detql=detql,detmf=tcmf,imergb=crgb, imemin=imemin,imemax=imemax,imesum=imesum,detbounds=detbounds) # sanity check for key in outkeys: if key not in outdict: print(key,'missing from returned dict') raw_input() return outdict def mf_contours(mf,mfdetcc,mfmapinfo,levels_ppmm,levels_rgba,level_max,mfmin=mfmin, mf_detids=[],mergedistm=50,scalef=2.0,detmeta={},doplot=False): 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 if levels_rgba.max()>1.0: levels_rgba = np.float32(levels_rgba)/255.0 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) np.set_printoptions(precision=9) if len(mf_detids)==0: mf_detids = np.unique(mfdetcc[mfdetcc!=0]) n_dets = len(mf_detids) contourmap = np.zeros([mf.shape[0],mf.shape[1]],dtype=np.float32) contourrgba = np.zeros([mf.shape[0],mf.shape[1],4],dtype=np.float32) mf_feats = [] if n_dets == 0: print('no detections, nothing to do') return contourmap, np.uint8(contourrgba*255), mf_feats if n_dets>1: print('processing %d detids: "%s"'%(n_dets,str(mf_detids))) elif mf_detids[0] not in (0,1): print('processing detid: "%s"'%str(mf_detids[0])) maxidx = np.zeros([3,len(mf_detids)],dtype=np.int32) levels_ordered = np.int32(sorted(levels_ppmm)) nlevel = len(levels_ordered) detmaxi = nlevel+1 selem_dilate = disk(3) selem_filter = disk(2) for detlab in mf_detids: detmask = mfdetcc==detlab (rmin,cmin),(rmax,cmax) = maskbbox(detmask,border=mfdist) mfsub = mf[rmin:rmax,cmin:cmax].copy() subr,subc = mfsub.shape[0],mfsub.shape[1] masksub = detmask[rmin:rmax,cmin:cmax] detmaxv = mfsub[masksub].max() if detmaxv < mfmin: print('detection id %d: no pixels above minppmm=%f'%(detlab,mfmin)) continue masksub = masksub*(mfsub>mfmin) ctrsub = contourmap[rmin:rmax,cmin:cmax] ctrrnz,ctrcnz = np.where(ctrsub!=0) # save r,c of maximum enhancement ccmaxrow,ccmaxcol = np.where(masksub & (mfsub==detmaxv)) ccmaxrow,ccmaxcol = ccmaxrow[0],ccmaxcol[0] if scalef != 1.0: rsr,rsc = int(round(subr*scalef)),int(round(subc*scalef)) mfsub = imresize(mfsub,[rsr,rsc],order=1) ctrsub = imresize(ctrsub,[rsr,rsc],order=0) masksub = imresize(mfsub,[rsr,rsc],order=0) ctridx = np.zeros_like(ctrsub,dtype=np.int32) ctrlonlat = [] #print('cclab',cclab,'npts',ccmask.sum(),'max',detmaxv) detlevels = levels_ordered[levels_ordered<=detmaxv] masksub = np.bool8(masksub!=0) # fill contour isolevels for li,ppmmi in enumerate(detlevels): detlevelmask = binary_fill_holes(masksub & (mfsub>=ppmmi)) detlevelmask = bwdilate(detlevelmask,selem=selem_dilate) detlevelmask = medianfilt(detlevelmask,selem_filter) detlevelr,detlevelc = np.where(detlevelmask!=0) if not detlevelmask.any(): print('no pixels found above level %d = %f ppmm'%(li,ppmmi)) continue ctrsub[detlevelr,detlevelc] = ppmmi ctridx[detlevelr,detlevelc] = li detlevelmask = np.float32(detlevelmask) 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 break else: print('no contours found at level %.3f'%(flevel)) print('contours for id=',detlab,'leveli=',li,'ppmmi=',ppmmi, 'num_contours=',len(cclevelctr), '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)) if doplot: pl.ioff() fig,ax = pl.subplots(1,3,sharex=True,sharey=True,num=1) ax[0].imshow(mfsub); ax[0].set_title('mfsub') ax[1].imshow(np.uint8(masksub)); ax[1].set_title('masksub') ax[2].imshow(ctridx,vmin=0,vmax=nlevel); ax[2].set_title('ctrsub') print(detlab,detmaxv,rmin,rmax,cmin,cmax,scalef) pl.show() # scale ctrsub back to original subimage size with interpolation to get a smooth image ccrme,cccme = int(ccmaxrow*scalef),int(ccmaxcol*scalef) ctrrgba = levels_rgba[ctridx].reshape([rsr,rsc,4]) if scalef != 1.0: ctrrgba = imresize(ctrrgba,[subr,subc,4],order=1,anti_alias=True) ctridx = imresize(ctridx,[subr,subc],order=0) ctridx = np.int32(np.round(ctridx)) ccrme,cccme = int(ccmaxrow),int(ccmaxcol) ctrrgba[ccrme,cccme,:] = level_max draw_crosshair=0 if draw_crosshair: chrad=4 ctrrgba[ccrme-chrad:ccrme+chrad+1,cccme,:] = level_max ctrrgba[ccrme,cccme-chrad:cccme+chrad+1,:] = level_max contourrgba[rmin:rmax,cmin:cmax,:] = ctrrgba # convert index values to ppmm, add maximum enhancement to contour map ctrlvl = levels_ppmm[ctridx] contourmap[rmin+ctrrnz,cmin+ctrcnz] = ctrlvl[ctrrnz,ctrcnz] contourmap[rmin+ccmaxrow,cmin+ccmaxcol] = detmaxv ccmaxlat,ccmaxlon = sl2latlon(rmin+ccmaxrow,cmin+ccmaxcol, mapinfo=mfmapinfo) ccprop = dict(levels=['%d'%l for l in detlevels], maxenhance=str((detmaxv,[ccmaxlat[0],ccmaxlon[0]]))) if detlab in detmeta: ccprop = detmeta[detlab] ccfeat = Feature(properties=ccprop,geometry=MultiPolygon(ctrlonlat)) mf_feats.append(ccfeat) if doplot: fig,ax = pl.subplots(1,3,sharex=True,sharey=True,num=1) ax[0].imshow(mfdetcc[rmin:rmax,cmin:cmax]); ax[0].set_title('mfdetcc') ax[1].imshow(contourmap[rmin:rmax,cmin:cmax]); ax[1].set_title('contourmap') ax[2].imshow(ctrrgba); ax[2].set_title('ctrrgba') pl.suptitle('detection id=%d'%detlab) pl.show() contourrgba = np.uint8(contourrgba*255) mf_feats = FeatureCollection(mf_feats) return contourmap, contourrgba, mf_feats def dtstr(): import datetime as dtime _now = dtime.datetime.now() return _now.strftime('%a %B %d, %Y %H:%M:%S') def process_mfimg(mfinfile,outdir,srcdata,colmap,args,NODATA=NODATA): mergedists = args.mergedists overwrite = args.clobber doplot = True #args.plot maxfetch = args.fetchmax nodetfilt = args.nodetfilt savefilt = True savefilt_only = args.savefilt_only or False colorf = args.colorfile mfmin = args.ppmmthr mfmax = args.ppmmmax skipdirs=False if skipdirs 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') return if not pathexists(outdir): mkdir(outdir) # default ppmm levels + colormap level_max = np.float32([ 0, 0, 0, 255])/255.0 # max enhancement = black by default levels_ppmm = np.float32([0,1000,2000,3000,4000]) levels_rgba = np.uint8([ [ 0, 0, 0, 0], # 0: < 0: transparent background [ 96, 96, 224, 255], # 0: darkblue [ 24, 224, 0, 255], # 2: green [200, 224, 0, 255], # 3: yellow [248, 0, 16, 255], # 4: red ]) levels_rgba = np.float32(levels_rgba)/255.0 #levels_hex = map(rgb2hex,levels_rgba[:,:3]/255.0) if pathexists(colorf): try: levels_ppmm, levels_rgba = parse_colors(colorf) if levels_ppmm[-1] == np.inf: level_max = levels_rgba[-1] levels_rgba = levels_rgba[:-1] levels_ppmm = levels_ppmm[:-1] except: print('Unable to parse color file "%s", using defaults'%colorf) pass lineids = srcdata[colmap['lid']] imgid = filename2flightid(mfinfile) imgdate = filename2flightdate(mfinfile) srcplumes = srcdata[lineids==imgid] mfdir,mffile = pathsplit(mfinfile) srcdirs = [pathjoin(outdir,'mergedist%d'%int(dist)) for dist in mergedists] srcdirs_exist = all([pathexists(srcdir) for srcdir in srcdirs]) statsf = pathjoin(outdir,mffile+'_ime.txt') if pathexists(statsf) and srcdirs_exist and not overwrite: print(statsf,'and all output directories already exist, nothing to do') print('(use --clobber argument to overwrite ALL existing output in %s)'%outdir) 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 if imgid.startswith('ang'): load_bands = [] # load everything rgb_bands = [0,1,2] elif imgid.startswith('f'): load_bands = [1,0] # load in reverse to get [gray,cmf] image rgb_bands = [1,1,1] # band 1 = cmf after reverse load mfdata = loadmaskedimage(mfinfile,load_bands=load_bands,masked_value=0, rgb_bands=rgb_bands) ch4mf = mfdata['image'] print('extrema(ch4mf): "%s"'%str((extrema(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' outdet_suf = out_suf+'_%d_%d'%(mfmin,mfmax) outkde_suf = out_suf+'_filt_det_kde_%d_%d'%(mfmin,mfmax) outcomp_suf = out_suf+'_filt_det_ccomp_%d_%d'%(mfmin,mfmax) out_base = splitext(pathsplit(mfinfile)[1])[0] det_outf = pathjoin(outdir,out_base+outdet_suf) kde_outf = pathjoin(outdir,out_base+outkde_suf) ccomp_outf = pathjoin(outdir,out_base+outcomp_suf) nothr_outf = pathjoin(outdir,out_base+out_suf) skip_kde = False areamin = args.minarea if not nodetfilt else 0 ch4det,detcc = filtdet(ch4mf,mfmask,mfmap,minarea=areamin, mfmin=mfmin,mfmax=mfmax,skip_kde=skip_kde, det_outf=det_outf,kde_outf=kde_outf, ccomp_outf=ccomp_outf) #detids = ccomp2detids(detcc,mf_ps) #point_ids = detids['point_ids'] #diffuse_ids = detids['diffuse_ids'] if savefilt: det_base = pathsplit(det_outf)[1] det_pngf = pathjoin(outdir,det_base+'.png') detrgb_pngf = pathjoin(outdir,det_base+'_rgb.png') # img_det = 32bit float mapped to 24-bit, 3-chan uint8, last channel dropped try: img_det = np.uint32(((2**24)-1)*ch4det).view(dtype=np.uint8) img_det = img_det.reshape([ch4det.shape[0],ch4det.shape[1],4]) img_det[...,3] = np.where(mfmask,0,255) imsave(det_pngf,np.uint8(img_det),verbose=True) #img_rgbql = rgbdet2ql(ch4rgb,ch4det,(detcc!=0)) #imsave(detrgb_pngf,np.uint8(img_rgbql),verbose=True) except: print("An error occurred saving %s"%det_pngf) pass if savefilt_only: return 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() lids = srcplumes[colmap['lid']].values 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 det_ids = {} for i,dist in enumerate(sorted(mergedists)): srcdir = srcdirs[i] if not pathexists(srcdir): print('creating product directory',srcdir) mkdir(srcdir) print('Merging adjacent ccomp ids for mergedist=%6.3f'%dist) dprev = mergedists[max(0,i-1)] if dprev in det_ids: # use previously merged ids computed using current distance threshold 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)))) if savefilt: ids_outf = pathjoin(srcdir,det_base+'_ids%d'%dist) array2img(ids_outf,ids,mapinfostr=mfmapstr,overwrite=True) det_ids[dist] = ids 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) ndists = len(mergedists) for i in range(nsrc): lid,sid,cid,pid = lids[i],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 = [NODATA]*ndists,[NODATA]*ndists detidi = [int(NODATA)]*ndists for disti,dist in enumerate(mergedists): idmask = det_ids[dist] print('Processing mergedist=',dist,'detections=',idmask.max()) imeout = compute_ime(mfdata,idmask,lat,lon,dist,maxfetch, mfmin=mfmin) 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' #ppmimgf = srcbase+'_cmf%d_%d.png'%(mfmin,mfmax) imeimgf = srcbase+'_ime.png' ctrimgf = srcbase+'_ctr.png' ctrgtif = ctrimgf.replace('.png','.tif') ctroutlimgf = srcbase+'_ctr_outl.png' figimgf = srcbase+'_fig.pdf' ctrf = srcbase+'_ctr.geojson' detrgb = imeout['detrgb'] #imsave(rgbimgf,detrgb,verbose=True) rgbgtif = rgbimgf.replace('.png','.tif') tile2geotiff(detrgb,detpos,rgbgtif,mfinfile) #detql = imeout['detql'] #qlimgf = srcbase+'_rgbql%d_%d.png'%(mfmin,mfmax) #imsave(qlimgf,detql,verbose=True) #qlgtif = qlimgf.replace('.png','.tif') #tile2geotiff(detql,detpos,qlgtif,mfinfile) detmf = imeout['detmf'] detmask = (detmf>=mfmin) #detql = rgbdet2ql(detrgb,detmf,detmask) #imsave(qlimgf,detql,verbose=True) if 0: if len(ppmql)!=0: print(ppmql), raw_input() 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 detdict = dict(Fetch=imeout['fetchm'], IME=imeout['imesum'], CCminppmm=mfmin, CCimevals=list(imeout['imevals']), CCrads=list(imeout['imerads']), CCcounts=list(imeout['cmfnums']), CCcmfvals=list(imeout['cmfvals'])) detmeta = {detid:detdict} mf_detids = np.int32([detid]) mf_ctr, mf_ctr_rgba, mf_feats = mf_contours(ch4mf,idmask,mfmap, levels_ppmm, levels_rgba, level_max, mf_detids=mf_detids, mfmin=mfmin, detmeta=detmeta) with open(ctrf,'w') as fid: print(str(mf_feats),file=fid) print('Saved contours to',ctrf) 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]) ctr_smooth = 0 #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]) ctr_rgba = extract_tile(mf_ctr_rgba,detpos,detdim,cval=0) #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))) # FIXME (BDB, 02/04/18): I believe this transpose is accurate, but need to verify # ctr_rgba = ctr_rgba.transpose(1,0,2) tile2geotiff(ctr_rgba,detpos,ctrgtif,mfinfile) #imsave(ctrimgf,ctr_rgba,verbose=True) #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[disti] = imeout['imesum'] fetchmi[disti] = imeout['fetchm'] detidi[disti] = imeout['detid'] ppmql = imeout['imergb'] imsave(imeimgf,ppmql,verbose=True) if doplot: fig,ax = pl.subplots(1,3,sharex=True,sharey=True,figsize=(12,5)) #ax[0].imshow(detql) ax[0].imshow(detrgb) ax[0].set_xlabel('lat,lon=%-6.3f,%-6.3f'%(plat,plon)) ax[1].imshow(ctr_rgba) ax[1].set_xlabel('IME$\Delta$=%4.2fm: [%-6.3f,%-6.3f]'%(imeout['radincm'], imeout['imemin'], imeout['imemax'])) ax[2].imshow(ppmql) ax[2].set_xlabel('PPM x m: [%-6.3f,%-6.3f]'%(mfmin,mfmax)) pl.suptitle('%s IME=%.3fkg, Fetch=%7.2fm (%.2f pixels)'%(cid,imeout['imesum'],imeout['fetchm'],imeout['fetchm']/mf_ps)) pl.setp([a.get_xticklabels() for a in fig.axes], visible=False) pl.savefig(figimgf) pl.clf() print('saved',figimgf) # TODO (BDB, 11/03/17): write outfiles for individual mergedists, # and combine afterwards 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) def find_cmf_file(ortdir,lid,calids,cmf_suffix): """ find_cmf_file(ortdir,lid,calids,cmf_suffix) Summary: summary Arguments: - ortdir: directory where orthorectified cmf images reside - lid: flightline id - calids: list of valid calibration ids - cmf_suffix: Output: - output """ cmf_regexp = pathjoin(ortdir,lid+cmf_suffix) mfmatch = glob(cmf_regexp) if len(mfmatch)==0: print('No CMF image found in directory',ortdir, 'matching line id:',lid) return None mfinfile = None foundcals = [] for mf in mfmatch: foundcals.append(filename2calid(mf)) for calid in calids: if not mfinfile and mf.endswith(calid+cmf_suffix.split('*')[1]): mfinfile = mf break if not mfinfile or not pathexists(mfinfile): print('No CMF image found in directory',ortdir, 'matching line id:',lid,'calids:',calids) print('Found the following calids for line', lid+':',foundcals) return None return mfinfile if __name__ == '__main__': import argparse parser = argparse.ArgumentParser(description="CMF IME Computation + Product Generation") parser.add_argument('-m','--mergedists', nargs='+', type=float, default=mergedists, help='List of one or more distances (in meters) to merge neighboring blobs (default=%s)'%mergediststr) #parser.add_argument('-d','--detfilt', type=str, # help='Path to filtered detection image (computed from mfinput by default)') #parser.add_argument('-p','--plot', action='store_true', # help='Save IME plots to outdir') parser.add_argument('-c','--clobber', action='store_true', help='Clobber (overwrite) existing outputs') 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('--minarea', type=int, help='Minimum area in pixels to remove small connected components in CMF (default=ang=%d, avcl=%d)'%(minarea_ang,minarea_avcl)) parser.add_argument('-p','--ppmmthr', type=float, default=mfmin, help='Minimum ppmm threshold (default=%.2f)'%(mfmin)) parser.add_argument('--ppmmmax', type=float, default=mfmax, help='Maximum ppmm threshold (default=%.2f)'%(mfmax)) parser.add_argument('--nodetfilt', action='store_true', help='Do not filter out small or weakly connected detections') parser.add_argument('-i','--inputfile', type=str, metavar='inputfile', help='CMF image (assumed listed in spreadsheet)') parser.add_argument('-l','--listfile', type=str, metavar='listfile', help='List of CMF images to process (subset of images in provided spreadsheet)') parser.add_argument('--calids', nargs='+', type=str, default=calids, help='List of one or more calibration IDs to consider (default=%s)'%str(calids)) parser.add_argument('--randomize', action='store_true', help='Randomize spreadsheet processing order (by source id)') parser.add_argument('--filter', type=str, metavar='filter', help='Compute IME for spreadsheet rows matching filters of the form "spreadsheet_column:filter_string (e.g., "flightline:ang201605")"') parser.add_argument('--savefilt_only', action='store_true', help='Save filtered detection images to outdir and exit (do not compute IME products)') parser.add_argument('--halt_errors', action='store_true', help='Stop processing flightlines if an error occurs (default=False)') parser.add_argument('--colorfile', type=str, metavar='colorfile', default=colorf, help='Color file (default="%s")'%colorf) 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 calids = args.calids inputfile = args.inputfile overwrite = args.clobber listfile = args.listfile xlsfile = args.xlsfile sheetname = args.sheetname srcfilter = args.filter randomize = args.randomize halt_errors = args.halt_errors srcdata = read_sources(xlsfile,sheetname,filterstr=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] ferr = None errorf = pathjoin(outdir,'error.log') nfiles = len(lineids) for i in range(nfiles): lid = lineids[i] outputdir = pathjoin(outdir,lid) if lid.startswith('ang'): cmf_dir = pathjoin(cmfdir_ang,'y'+lid[5:7],'cmf','ort') cmf_suffix = cmfsuf_ang args.minarea = args.minarea or minarea_ang elif lid.startswith('f'): cmf_dir = cmfdir_avcl cmf_suffix = cmfsuf_avcl args.minarea = args.minarea or minarea_avcl cmffile = find_cmf_file(cmf_dir,lid,calids,cmf_suffix) if not cmffile: print('CMF file for line "%s" not found'%lid) continue tstamp = dtstr() try: print('Processing cmffile=',cmffile, '(%d of %d)'%(i+1,nfiles)) print('start time=',tstamp) process_mfimg(cmffile,outputdir,srcdata,colmap,args) print('-'*65) print('') except KeyboardInterrupt: print('Recieved keyboard interrupt, exiting') break except: import traceback if ferr is None: ferr = open(errorf,'w') else: print('-'*65,file=ferr) print('[%s] %s errors:'%(tstamp,lid),file=ferr) exc_msg = '[%s] An unexpected exception of type %s occurred%s: %s' exc_type, value, tb = sys.exc_info() exc_file = tb.tb_frame.f_locals.get('filename',None) msg_file = ' in file %s'%exc_file if exc_file else '' # log to stdout err_msg = exc_msg % (tstamp,str(exc_type),msg_file,str(value)) print(err_msg) print('[%s] Traceback:'%tstamp) traceback.print_tb(tb,None) # log to error.log print(err_msg,file=ferr) print('[%s] Traceback:'%tstamp,file=ferr) traceback.print_tb(tb,None,ferr) print('[%s] Errors logged to %s'%(tstamp,errorf)) if halt_errors: break if ferr: ferr.close()