Python util.mk_dir函数代码示例

        if binPr==0.0:
          a2cltmp  = ma.masked_where(a2ra>binPr,  a2cl)
        else:
          a2cltmp  = ma.masked_where(a2ra>=binPr, a2cl)

        for cltype in lcltype:
          clid    = dclid[cltype]
          a2tmp = ma.masked_where(a2cltmp !=clid, a2pr)
    
          da3sum[binPr][cltype] = da3sum[binPr][cltype] + a2tmp.filled(0.0)
          da3num[binPr][cltype] = da3num[binPr][cltype] + ma.masked_where(a2tmp.mask, a2oneint).filled(0)
  
    # Save Monthly output
    #baseDir = "/tank/utsumi/PMM/WNP.261x265"
    baseDir = "/home/utsumi/mnt/well.share/PMM/WNP.261x265"
    oDir    = baseDir + "/ByRA.CL.Pr.%s"%(dattype)
    util.mk_dir(oDir)
  
    # Bin file
    binPath = oDir + "/CloudType.txt"
    sout    = "\n".join(["%d:%d"%(cltype, dclid[cltype]) for cltype in lcltype]).strip()
    f=open(binPath,"w"); f.write(sout); f.close()
  
    for binPr in lbinPr:
      sumPath   = oDir + "/sum.P%05.1f.%04d.%02d.%dx%dx%d"%(binPr, Year,Mon, ncltype,ny,nx)
      numPath   = oDir + "/num.P%05.1f.%04d.%02d.%dx%dx%d"%(binPr, Year,Mon, ncltype,ny,nx)
  
      da3sum[binPr].tofile(sumPath)
      da3num[binPr].tofile(numPath)  
      print sumPath

            plt.loglog(agvTmp, aesurfTmp, 'o')
            
            plt.ylim([0.01,100])
            plt.xlim([0.01,100])
            plt.savefig(figPath)
            print figPath
            plt.clf()


#--- save data ----
for prtype in lprtype:
    if calc != True:
        continue

    outDir = '/home/utsumi/mnt/wellshare/GPMGV/dt-eSurf-%s/dist.%.1fkm.ndom.%02d.%04d.%02d-%04d.%02d'%(prdName, thdist, len(ldomain), iYM[0], iYM[1], eYM[0], eYM[1])
    util.mk_dir(outDir)

    for dattype in ldattype:  
        datPath= outDir + '/dt-eSurf.%s.%s.%s.npy'%(prdName, dattype, prtype)

        if dattype=='cc':
            a1dat = da1cc[prtype]
        elif dattype=='rmse':
            a1dat = da1rmse[prtype]
        elif dattype=='bias':
            a1dat = da1bias[prtype]
        elif dattype=='brat':
            a1dat = da1brat[prtype]
        elif dattype=='num':
            a1dat = da1num[prtype]
        elif dattype=='rain':

        dpr[lndsea,icl].extend(a2prtmp.compressed())
        dku[lndsea,icl].extend(a2kutmp.compressed())

    #-- Save --------- 
    lkeys = [[lndsea,icl] for lndsea in llndsea
                          for icl     in lcltype]
    for [lndsea,icl] in lkeys:
      #baseDir  = "/home/utsumi/mnt/wellshare/PMM/WNP.261x265"
      baseDir  = ibaseDir
      sDir     = baseDir + "/VsKuPR.CL.%s/%04d"%(dattype,Year)

      prPath   = sDir + "/%s.%04d.%02d.%s.%s.bn"%(dattype,Year,Mon,lndsea, dclShortName[icl])
      kuPath   = sDir + "/KuPR.%04d.%02d.%s.%s.bn"%(Year,Mon,lndsea, dclShortName[icl])

      util.mk_dir(sDir)
      array(dpr[lndsea,icl], float32).tofile(prPath)
      array(dku[lndsea,icl], float32).tofile(kuPath)

      print prPath

    # Combine lnd, sea, cst
    for icl in lcltype:
      baseDir = ibaseDir
      sDir  = baseDir + "/VsKuPR.CL.%s/%04d"%(dattype,Year)

      lstype = [dattype, "KuPR"]

      for stype in lstype:
        lndPath = sDir + "/%s.%04d.%02d.%s.%s.bn"%(stype,Year,Mon,"lnd", dclShortName[icl])
        seaPath = sDir + "/%s.%04d.%02d.%s.%s.bn"%(stype,Year,Mon,"sea", dclShortName[icl])

    urlon = float( Decimal(lon1).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP) )

    a1lat_map = arange(lllat,urlat+res*0.1, res)
    a1lon_map = arange(lllon,urlon+res*0.1, res)

    for YM in lYM:
        Year   = YM[0]
        Mon    = YM[1]

        if (domain,Year,Mon) not in dgName.keys():
            print 'no obs',domain,Year,Mon
            continue
        #-- distance from nearest gauge -
        a2mindist = ret_a2mindist(domain, Year, Mon)
        gvmapDir  = gvmapbaseDir + '/%s/%04d%02d'%(domain,Year,Mon)
        util.mk_dir(gvmapDir)
        distPath  = gvmapDir + '/map.mindist.%s.%04d%02d.npy'%(domain,Year,Mon)
        np.save(distPath, a2mindist)
        print 'make minimum distance file'
        print distPath

        #-- check satellite overpass time
        sateDir = satebaseDir + '/%s/%04d%02d'%(domain, Year, Mon) 
        if not os.path.exists(sateDir):
            print 'No directory', sateDir
            continue
        ssearch = sateDir + '/prcp.*.npy'
        lsatePath = glob.glob(ssearch)
        lsatePath = sorted(lsatePath)

#simbaseDir  = '/home/utsumi/work/temp/outflw'
#loclistDir  = '/home/utsumi/work/temp/global_15min'
#loclistPath = loclistDir + '/grdc_loc.txt'
#outDir      = '/home/utsumi/work/temp/timeseries'
#figDir      = '/home/utsumi/work/temp/fig'

grdcDir     = '/work/data2/GRDC/dat_day'
simbaseDir  = '/work/hk01/utsumi/PMM/hydro/outflw.trip'
#loclistPath = '/work/hk01/utsumi/PMM/hydro/list/grdc_loc_rev_1998_2009.txt'
outDir      = '/work/hk01/utsumi/PMM/hydro/timeseries.trip'
figDir      = '/work/hk01/utsumi/PMM/hydro/fig'




util.mk_dir(outDir)
util.mk_dir(figDir)
#*** Functions *********************
def load_sim(prcp, stnid, iYear, eYear):
    nYear = eYear - iYear + 1
    lYear = range(iYear,eYear+1)
    a1out = deque([])
    for Year in lYear:
        simDir  = simbaseDir + '/%s/%04d'%(prcp, Year)
        srcPath = simDir + '/flwout.%07d.npy'%(stnid)
        a1tmp   = np.load(srcPath)
        a1out.extend(a1tmp) 

    return array(a1out)

def get_timerange_grdc(srcPath):

  return [s.strip() for s in lines]

#******************************

ltag   = ret_ltag(ntag)
for Year,Mon in lYM:
  eDay   = calendar.monthrange(Year,Mon)[1]
  iDTime = datetime(Year,Mon,1,0)
  eDTime = datetime(Year,Mon,eDay,23)
  dDTime = timedelta(hours=24)
  lDTime = util.ret_lDTime(iDTime, eDTime, dDTime)

  a3tagprcp = zeros([ntag, ny, nx], float32)
  for DTime in lDTime:
    print DTime
    a2prcp  = load_Prcp(DTime, PRJ)

    frcPath = CH08.path_TagH08(DTime, model="JRA55", wnflag="nn", nclass=ntag, res="one", tstep="day")[1]
    a3frc   = fromfile(frcPath, float32).reshape(ntag, ny, nx)
    for i in range(ntag):
      a3tagprcp[i] = a3tagprcp[i] + a2prcp * a3frc[i] 

  a3tagprcp = a3tagprcp / len(lDTime)   # average
  #-- save -----
  rootDir, srcDir, srcPath = h08.ret_pathTagVarMon(var="PRCP", ntag=ntag, Year=Year, Mon=Mon) 

  util.mk_dir(srcDir)
  a3tagprcp.tofile(srcPath)
  print srcPath 

grdc = GRDC.GRDC()

iYear = 2000
eYear = 2009

lprcp  =['GPCC','GPCC','GPCC']
lmetName=['CC','NSE','RBias']
#calc  = True
calc  = False

grdcDir     = '/work/data2/GRDC/dat_day'
simbaseDir  = '/work/hk01/utsumi/PMM/hydro/outflw.trip'
outDir      = '/work/hk01/utsumi/PMM/hydro/metric.trip'
csvDir      = '/work/hk01/utsumi/PMM/hydro' + '/csv.trip'
util.mk_dir(outDir)
util.mk_dir(csvDir)
#*** Functions *********************
def load_sim(prcp, stnid, iYear, eYear):
    nYear = eYear - iYear + 1
    lYear = range(iYear,eYear+1)
    a1out = deque([])
    for Year in lYear:
        simDir  = simbaseDir + '/%s/%04d'%(prcp, Year)
        srcPath = simDir + '/flwout.%07d.npy'%(stnid)
        a1tmp   = np.load(srcPath)
        a1out.extend(a1tmp) 
    return array(a1out)

def get_timerange_grdc(srcPath):
    f=open(srcPath,'r'); lines=f.readlines(); f.close()

  a2ave1 = ma.masked_where(a2num==0.0, a2sum1)/a2num
  a2ave2 = ma.masked_where(a2num==0.0, a2sum2)/a2num

  #****** Figure ***************
  figplot  = plt.figure(figsize=(5,5))
  axplot   = figplot.add_axes([0.13, 0.13, 0.77, 0.77])
  axplot.plot( a2ave1, a2ave2, "o", color="k", markersize=2)

  #-- axis -----
  axmax = 15.
  axplot.set_xlim(0.0, axmax)
  axplot.set_ylim(0.0, axmax)
 
  #-- axis label --
  axplot.set_xlabel("%s [mm/hour]"%(prtype1), fontsize=18)
  axplot.set_ylabel("%s [mm/hour]"%(prtype2), fontsize=18)

  #-- title -------
  Corr   = corrcoef(a2ave1.flatten(), a2ave2.flatten())[0][1]
  stitle = "X=%s Y=%s tag:%s"%(prtype1, prtype2,tag)
  stitle = stitle + "\n%04d-%04d M%02d-%02d R=%5.2f"%(iYear,eYear, lMon[0], lMon[-1], Corr)
  axplot.set_title(stitle, fontsize=15)

  #-- save -------
  figDir  = baseDir + "/pict"
  figPath = figDir  + "/scatter.%04d-%04d.M%02d-%02d.%s.png"%(iYear,eYear,lMon[0],lMon[-1],tag)
  util.mk_dir(figDir)
  plt.savefig(figPath)
  print figPath

    continue


  nl, nw        = gpmobt.data.shape
  a2sum, a2num  = pmm_fsub.obt2wnpac261x265(gpmobt.data.T, gpmobt.lon.T, gpmobt.lat.T, nw, nl)

  a2sum   = a2sum.T
  a2num   = a2num.T
  a2pr    = ma.masked_invalid(a2sum/a2num)
  #a2pr    = a2pr / (60.*60.)   # mm/h --> mm/s
  a2pr    = a2pr.filled(miss_out)
  #--- save ---------
  Year    = DTime.year
  Mon     = DTime.month
  Day     = DTime.day
  Hour    = DTime.hour
  Minute  = DTime.minute
  baseDir = "/tank/utsumi/PMM/WNP.261x265/%s/%s/%s"%(prj,prdLv,prdVer)
  dataDir = baseDir + "/%04d/%02d"%(Year,Mon)
  dataPath= dataDir + "/pr.%02d.%02d.%02d.%dx%d"%(Day,Hour,Minute,ny,nx)

  util.mk_dir(dataDir)
  a2pr.astype(float32).tofile(dataPath)
  print "OUTPUT"
  print dataPath
    
  

    asx = zeros([nh,nt])
    asy = zeros([nh,nt])
    asxx= zeros([nh,nt])
    asyy= zeros([nh,nt])
    asxy= zeros([nh,nt])
    anum= zeros([nh,nt])

    for domain in ldomain:
        for YM in lYM:
            Year, Mon = YM
            if (domain,Year,Mon) not in dgName.keys():
                print 'no obs',domain,Year,Mon
                continue

            datDir = '/home/utsumi/mnt/wellshare/GPMGV/dt-lev-PR/%s'%(domain)
            util.mk_dir(datDir)

            satePath    = datDir + '/ssate.%04d%02d.%s.npy'%(Year,Mon,prtype)
            sate2Path   = datDir + '/ssate2.%04d%02d.%s.npy'%(Year,Mon,prtype)
            gvPath      = datDir + '/sgv.%04d%02d.%s.npy'%(Year,Mon,prtype)
            gv2Path     = datDir + '/sgv2.%04d%02d.%s.npy'%(Year,Mon,prtype)
            sategvPath  = datDir + '/sategv.%04d%02d.%s.npy'%(Year,Mon,prtype)
            numPath     = datDir + '/num.%04d%02d.%s.npy'%(Year,Mon,prtype)

            asxTmp  = np.load(gvPath    )
            asxxTmp = np.load(gv2Path   )
            asyTmp  = np.load(satePath  )
            asyyTmp = np.load(sate2Path )
            asxyTmp = np.load(sategvPath)
            anumTmp = np.load(numPath   )

#.........这里部分代码省略.........
      ytick = None

    elif figtype=="cnt.focus":
      xtick = arange(0.0, 6.0+0.01, 1.0)
      sxtick = None 
      ytick = None
  
    elif figtype =="Weak.pdf":
      xtick = arange(0, 3+0.01, 0.5)
      sxtick = None 
  
    elif figtype =="Heavy.pdf":
      if icltype in [1,2]:
        xtick = arange(0, 0.05+0.001, 0.01)
        sxtick = None 
      else:
        ytick = [0, 1e-5, 1e-4, 1e-3]
  
      axplot.set_xlim(10, lbinPr[-2])
      xtick = arange(10, lbinPr[-2]+1, 10) 
      sxtick = None 
    elif figtype =="log.Weak.pdf":
      plt.yscale("log")
      if icltype in [1,2]:
        ytick = [1.e-3, 1.e-2, 1.e-1, 1.e+0]
        xtick = arange(0, 3+0.01, 0.5)
        sxtick = None 
      else:
        ytick = [1.e-6, 1.e-5, 1.e-4, 1.e-3, 1.e-2, 1.e-1, 1.e+0, 1.e+1]
        xtick = arange(0, 3+0.01, 0.5)
        sxtick = None 
  
    elif figtype =="log.Heavy.pdf":
      plt.yscale("log")
      if icltype in [1,2]:
        ytick = [1.e-5, 1.e-4, 1.e-3, 1.e-2, 1.e-1]
      else:
        axplot.set_ylim(1.e-8, 1.e-3)
        ytick = [1.e-8, 1.e-7, 1.e-6, 1.e-5, 1.e-4, 1.e-3]
      xtick = arange(10, lbinPr[-2]+1, 10) 
      sxtick = None 
  else:
    print "check figtype", figtype
    sys.exit()

  # Axis-tick labels
  if xtick !=None:
    axplot.set_xlim([xtick[0],xtick[-1]])
    if sxtick !=None:
      plt.xticks( xtick, sxtick)
    else:
      plt.xticks( xtick, xtick)
  if ytick !=None:
    axplot.set_ylim([ytick[0],ytick[-1]])
    plt.yticks( ytick, ytick)

  for tick in axplot.xaxis.get_major_ticks():
    tick.label.set_fontsize(19)
  for tick in axplot.yaxis.get_major_ticks():
    tick.label.set_fontsize(16)

  # Add title
  stitle = "%s %04d/%02d-%04d/%02d CL=%s"\
             %(lndsea,iYM[0],iYM[1],eYM[0],eYM[1],dclName[icltype])

  if   expr=="old":
    stitle = stitle + " old"
  elif expr=="sht":
    stitle = stitle + " sht"


  if dommask !=None:
    stitle = stitle + " %sdom"%(dommask)

  plt.title(stitle, fontsize=10) 

  # Save figure
  if expr  =="std":
      sDir  = ibaseDir + "/pict"
  elif expr=="old":
      sDir  = ibaseDir + "/pict.old"
  elif expr=="sht":
      sDir  = ibaseDir + "/pict.sht"
  util.mk_dir(sDir)
  if dommask==None:
    sPath = sDir + "/%s.%s.%s.png"%(figtype, lndsea, dclShortName[icltype]) 
  else:
    sPath = sDir + "/%sdom.%s.%s.%s.png"%(dommask, figtype, lndsea, dclShortName[icltype]) 
  plt.savefig(sPath)
  print sPath

  plt.close()
  # Legend file
  legPath = sDir + "/legend.pdfs.png"
  figleg  = plt.figure(figsize=(3,3))
  lines = [line[0] for line in lines]  # 2D list to 1D list
  figleg.legend(lines, ret_ldatname(ldattype), fontsize=15)
  figleg.savefig(legPath)
 
  plt.close()

BBox   = ra.BBox
LatOrg = ra.Lat
LonOrg = ra.Lon
LatUp  = arange(BBox[0][0]+resUp/2,
                BBox[1][0]+resUp/10,
                resUp, dtype='float64')

LonUp  = arange(BBox[0][1]+resUp/2,
                BBox[1][1]+resUp/10,
                resUp, dtype='float64')

nyUp, nxUp = len(LatUp), len(LonUp)
baseDir  = "/tank/utsumi/data/RadarAMeDAS/ra_%s"%(resUp)
#--- save meta data -----------
util.mk_dir(baseDir)
sLat = util.array2csv(LatUp)
sname= baseDir+"/lat.csv"
f=open(sname,"w");f.write(sLat); f.close()

sLon = util.array2csv(LonUp)
sname= baseDir+"/lon.csv"
f=open(sname,"w");f.write(sLon); f.close()
print sname

#------------------------------
def ret_oPath(DTime, resUp):
  Year, Mon, Day, Hour, Min =\
    DTime.year, DTime.month, DTime.day, DTime.hour, DTime.minute
  oDir     = os.path.join(baseDir, "%04d%02d"%(Year,Mon))
  oPath    = os.path.join(oDir, "RadarAmedas.%04d%02d%02d%02d%02d00.%dx%d"%(Year,Mon,Day,Hour,Min,nyUp, nxUp))

#.........这里部分代码省略.........
      if icltype==1:
        axplot.set_ylim(1.e-3, 1.e+0)
        axplot.set_xlim(0.0, 3.0)
      else:
        axplot.set_ylim(1.e-6, 1.e+1)
        axplot.set_xlim(0.0, 3.0)
  
    elif figtype =="log.Heavy.pdf":
      plt.yscale("log")
      if icltype==1:
        axplot.set_ylim(1.e-5, 1.e-1)
      else:
        axplot.set_ylim(1.e-8, 1.e-3)
      axplot.set_xlim(10.0, lbinPr[-2])
    

  elif clVer[:5] == "MyWNP":
    if figtype=="pdf":
      axplot.set_ylim(0.0, 0.3)
      axplot.set_xlim(0.0, 20.0)
 
    elif figtype == "pdf.focus":
      axplot.set_ylim(0.0, 0.05)
      axplot.set_xlim(0.0, 6.0)
 
    elif figtype=="cnt":
        axplot.set_xlim(0.0, 20.0)

    elif figtype=="cnt.focus":
        axplot.set_xlim(0.0, 6.0)
  
    elif figtype =="Weak.pdf":
      axplot.set_xlim(0.0, 3.0)
  
    elif figtype =="Heavy.pdf":
      if icltype in [1,2]:
        axplot.set_ylim(0.0, 0.05)
      else:
        axplot.set_ylim(0.0, 1e-3)
  
      axplot.set_xlim(10, lbinPr[-2])
  
    elif figtype =="log.Weak.pdf":
      plt.yscale("log")
      if icltype in [1,2]:
        axplot.set_ylim(1.e-3, 1.e+0)
        axplot.set_xlim(0.0, 3.0)
      else:
        axplot.set_ylim(1.e-6, 1.e+1)
        axplot.set_xlim(0.0, 3.0)
  
    elif figtype =="log.Heavy.pdf":
      plt.yscale("log")
      if icltype in [1,2]:
        axplot.set_ylim(1.e-5, 1.e-1)
      else:
        axplot.set_ylim(1.e-8, 1.e-3)
      axplot.set_xlim(10.0, lbinPr[-2])
     
  else:
    print "check figtype", figtype
    sys.exit()

  # Axis-tick labels
  for tick in axplot.xaxis.get_major_ticks():
    tick.label.set_fontsize(20)
  for tick in axplot.yaxis.get_major_ticks():
    tick.label.set_fontsize(17)

  # Add title
  stitle = "%04d/%02d-%04d/%02d CL=%s"\
             %(iYM[0],iYM[1],eYM[0],eYM[1],dclName[icltype])

  if dommask !=None:
    stitle = stitle + " %sdom"%(dommask)

  plt.title(stitle, fontsize=10) 

  # Save figure
  #sDir  = "/tank/utsumi/PMM/WNP.261x265/pict"
  #sDir  = "/home/utsumi/mnt/well.share/PMM/WNP.261x265/pict"
  #sDir  = "/home/utsumi/mnt/well.share/PMM/WNP.261x265/pict.MyCL"
  sDir  = ibaseDir + "/pict"
  util.mk_dir(sDir)
  if dommask==None:
    sPath = sDir + "/%s.%s.png"%(figtype, dclShortName[icltype]) 
  else:
    sPath = sDir + "/%sdom.%s.%s.png"%(dommask, figtype, dclShortName[icltype]) 
  plt.savefig(sPath)
  print sPath

  plt.close()
  # Legend file
  legPath = sDir + "/legend.pdfs.png"
  figleg  = plt.figure(figsize=(2,3))
  lines = [line[0] for line in lines]  # 2D list to 1D list
  figleg.legend(lines, ldattype)
  figleg.savefig(legPath)
 
  plt.close()