Changeset 1070 for trunk/UTIL/PYTHON

Timestamp:

Oct 14, 2013, 5:36:36 PM (11 years ago)

Author:

aslmd

Message:

PYTHON. updated analysis tools for dust devils.

Location:

trunk/UTIL/PYTHON/powerlaw

Files:

• ddhist.py (modified) (4 diffs)
• ddstat.py (modified) (1 diff)
• dropsize.py (modified) (1 diff)
• finddd.py (modified) (5 diffs)
• performdd.py (modified) (4 diffs)

trunk/UTIL/PYTHON/powerlaw/ddhist.py

@@ -5 +5 @@
 import scipy.optimize as sciopt
 
-## SAVED
-namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = False ; typefit=1
-namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = True ; typefit=2
-namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = False ; typefit=1
-namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = True ; typefit=1
-namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt" ; drop = False ; typefit=1 #bof. tous les 10 est mieux.
+### SAVED
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = False ; typefit=1
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = True ; typefit=2
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = False ; typefit=1
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = True ; typefit=1
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt" ; drop = False ; typefit=1 #bof. tous les 10 est mieux.
+#
+###########################################################
+##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_1.txt"
+##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_1.txt"
+##namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_1.txt"
+##namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc1.txt"
+#namefile = "/home/aymeric/Big_Data/LES_dd/sav/press_ustm_exomars.ncm2_1.txt"
+##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_1.txt"
 
-##########################################################
-#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_1.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_1.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_1.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt"
-namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc1.txt"
-namefile = "/home/aymeric/Big_Data/LES_dd/sav/press_ustm_exomars.ncm2_1.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_1.txt"
+
+
+
+case = "188324p"
+case = "191798"
+case = "160564p"
+case = "156487"
+case = "2007p"
+case = "13526p"
+case = "172097"
+namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_1.txt"
+
+
+
 drop = False ; typefit=1
 #drop = True ; typefit=1
@@ -25 +40 @@
 ##########################################################
 ## define bins (we expect fit do not depend too much on this -- to be checked)
-##nbins = 7 ; www = 3.
+nbins = 7 ; www = 3.
 #nbins = 10 ; www = 2.5
-nbins = 15 ; www = 2. 
+#nbins = 15 ; www = 2. 
 ##nbins = 20 ; www = 1.75
 ##nbins = 30 ; www = 1.5
 ##nbins = 50 ; www = 1.2
+##########################################################
+limrest = 4. # restrict limit (multiple of dx)
+#limrest = 2.
+#limrest = 3.
+#limrest = 0.
 ##########################################################
 
@@ -56 +76 @@
 # restrictions
 restrict = (s > 0) # initialization (True everywhere)
+restrict = restrict*(s >= limrest*dx) # remove lowest sizes (detection limit) 
 #restrict = restrict*(np.abs(i-j) <= 6.*dx) # condition sur i,j (width,height)
 #restrict = restrict*(d > 0.9) # limit on drop for size (casse la power law? seulement si keep smaller devils)
 #restrict = restrict*(d > 0.5)
-#restrict = restrict*(s >= 3.*dx) # remove lowest sizes (detection limit) 
 var = var[restrict]
 
@@ -76 +96 @@
 print "min %5.2e // max %5.2e" % (np.min(var),np.max(var))
 for iii in range(len(zebins)-1):
-    print "%5.2e in [%5.2e %5.2e]" % (yeah[0][iii],zebins[iii],zebins[iii+1])
+    print "%5.2e in [%5.0f %5.0f] %5.0f" % (yeah[0][iii],zebins[iii],zebins[iii+1],np.abs(zebins[iii+1]-zebins[iii]))
 
 # fitting function to superimpose

trunk/UTIL/PYTHON/powerlaw/ddstat.py

@@ -12 +12 @@
 #namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_2.txt"
 ##########################################################
+
+
+ 
+
+#case = "188324p"
+#case = "191798"
+#case = "160564p"
+#case = "156487"
+#case = "2007p"
+#case = "13526p"
+#case = "172097"
+
+
+namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_2.txt"
+
+
+
+
 
 # load data

trunk/UTIL/PYTHON/powerlaw/dropsize.py

@@ -13 +13 @@
 ##########################################################
 
-# read data np.loadtxt?
-## import pylab as plb data = plb.loadtxt('data.dat') x = data[:,0] y= data[:,1]
-f = open(namefile,'r')
-t = [] ; s = [] ; d = [] ; i = [] ; j = []
-for line in f:
-    ct,cs,cd,ci,cj = line.strip().split(';')
-    t.append(float(ct))
-    s.append(float(cs)) ; d.append(float(cd))
-    i.append(int(ci)) ; j.append(int(cj))
-f.close()
-# ensure numpy array
-t = np.array(t) ; s = np.array(s) ; d = np.array(d) ; i = np.array(i) ; j = np.array(j)
+#case = "188324p"
+#case = "191798"
+#case = "160564p"
+#case = "156487"
+#case = "2007p"
+#case = "13526p"
+#case = "172097"
 
-drop = ppplot.plot1d() ; drop.field = d ; drop.absc = s
-drop.lstyle = '' ; drop.marker = '.' ; drop.color = 'r' ; drop.fmt = "%.1f"
+
+namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_1.txt"
+
+# load data
+data = np.loadtxt(namefile,delimiter=";")
+t = data[:,0] ; s = data[:,1] ; d = data[:,2]
+i = data[:,3] ; j = data[:,4] # int?
+
+drop = ppplot.plot1d() ; drop.f = d ; drop.x = s
+drop.linestyle = '' ; drop.marker = '.' ; drop.color = 'r' ; drop.fmt = "%.1f"
 drop.xlabel = "Vortex size (m)" ; drop.ylabel = "Pressure drop (Pa)"
-drop.make() ; ppplot.save()
+drop.makeshow()
 

trunk/UTIL/PYTHON/powerlaw/finddd.py

@@ -1 +1 @@
 #! /usr/bin/env python
-from ppclass import pp
+from ppclass import pp,ncattr
 from ppplot import plot2d,save,writeascii
 import numpy as np
@@ -19 +19 @@
 ###############################################################################
 def finddd(filefile,\
-           timelist,\
-           dx=50.,\
+           timelist=None,\
            dt_out=50.,\
            lt_start=8.,\
-           halolim=30.,\
+           halolim=None,\
            method=1,\
            plotplot=False,\
@@ -52 +51 @@
         myfile2 = open(filefile+'m'+str(method)+'_'+'2.txt', 'w')
     ###############################################################################
-    
+
+    ## get the resolution within the file
+    dx = ncattr(filefile,'DX') ; print "resolution in meters is: ",dx
+    ## if no halolim is given, guess it from resolution
+    if halolim is None:
+        extentlim = 2000. # the putative maximum extent of a vortex in m
+        halolim = extentlim / dx
+        print "maximum halo size is: ",halolim
+
     ## mean and std calculations
     ## -- std is used in both methods for limits
@@ -72 +79 @@
     print "**************************************************************"
     
+    # if no timelist is given, take them all
+    if timelist is None:
+        sizet = psfc.shape[0]
+        print "treat all time values: ",sizet
+        timelist = range(0,sizet-1,1)
+
     ## LOOP ON TIME
     for time in timelist:
@@ -80 +93 @@
      #ustm = pp(file=filefile,var="USTM",t=time).getf()
 
-     ## MAIN ANALYSIS. LOOP ON FAC.
-     for fac in faclist:
-     #fac = 3.75
-     #for method in [2,1]:
+     ## MAIN ANALYSIS. LOOP ON FAC. OR METHOD.
+     #for fac in faclist:
+     fac = 3.75
+     for method in [2,1]:
   
       ## initialize arrays

trunk/UTIL/PYTHON/powerlaw/performdd.py

@@ -1 +1 @@
 #! /usr/bin/env python
 from finddd import finddd
+
+finddd("/planeto/aslmd/LESdata/case_A_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
+
+exit()
+
+# LES_DUST_DEVIL_OLDPHYSok
+finddd("/planeto/aslmd/LESdata/188324p.nc",dt_out=50.,lt_start=9.)
+# LES
+finddd("/planeto/aslmd/LESdata/160564p.nc",dt_out=50.,lt_start=9.)
+finddd("/planeto/aslmd/LESdata/156487.nc",dt_out=50.,lt_start=8.)
+# LES_INSIGHT
+finddd("/planeto/aslmd/LESdata/2007p.nc",dt_out=10.,lt_start=6.)
+finddd("/planeto/aslmd/LESdata/13526p.nc",dt_out=10.,lt_start=6.)
+# EXOMARSshear
+finddd("/planeto/aslmd/LESdata/172097.nc",dt_out=50.,lt_start=13.)
+# QJ cases
+finddd("/planeto/aslmd/LESdata/case_A_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
+finddd("/planeto/aslmd/LESdata/case_B_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
+finddd("/planeto/aslmd/LESdata/case_C_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
+finddd("/planeto/aslmd/LESdata/case_I_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
+finddd("/planeto/aslmd/LESdata/case_HIGH_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
+
+#####################################################################
+#####################################################################
+#####################################################################
+
 
 ####################
@@ -8 +34 @@
 ####################
 ####################
-
 
 #finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",range(0,591,10),method=2)
@@ -22 +47 @@
 ####################
 ####################
-finddd("/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc",range(0,54,1),dx=12.,halolim=100,method=2)
+#finddd("/planeto/aslmd/LESdata/press_ustm_exomars.nc",range(0,54,1),halolim=100,method=2)
 #timelist = [28.,48.,49.]
 #plotplot = True
@@ -32 +57 @@
 ####################
 #finddd("/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.nc",\
-#           range(0,336,1),dx=10.,lt_start=9.,halolim=150,method=2)
+#           range(0,336,1),lt_start=9.,halolim=150,method=2)
 ###halolim 50 idem
 ####################