Changeset 349 for trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples

Timestamp:

Nov 5, 2011, 6:54:37 PM (13 years ago)

Author:

aslmd

Message:

PYTHON GRAPHICS: A REFERENCE VERSION FOR BOTH GCM AND MESOSCALE USE!

Merged recent changes by TN for better GCM handling + 1D/2D plots with AS mesoscale version [winds.py]

Now

• there is a common function named planetoplot in planetoplot.py which supports both GCM and mesoscale plots

  this will be a common reference for future development
  

• mesoscale users might want to use the script meso.py
• GCM users might want to use the script gcm.py

Todo

• a few capabilities that were in winds.py are broken in meso.py, this will be fixed soon. the 'old' reference script winds.py is still here anyway!
• there is still a bug with 'ortho' projection with GCM fields. still searching...

Misc

• A few modifications were done in myplot.py and mymath.py
• Added a simple stupid little.py script illustrating very simple use of Python graphics
• Added the script doing statistical analysis for dust devils

File:

• trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py (modified) (8 diffs)

trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py

@@ -27 +27 @@
    
     ### LOAD NETCDF DATA
-    filename = "psfc.nc"
+    #filename = "/home/aymeric/Big_Data/psfc.nc"
     nc = Dataset(filename) 
     psfc = nc.variables["PSFC"]
@@ -36 +36 @@
     allsizesx = []
     allsizesy = []
+    depression = []
     stride = 1 #5
     for i in range(0,shape[0],stride):
@@ -52 +53 @@
         fac = 3.5
         #fac = 2.5
+        fac = 3.
         lim = ave - fac*std
         where = np.where(psfc2d < lim)
         ############### END CRITERION
+
+        depression = np.append(depression,np.ravel(psfc2d[where])-ave)
    
         lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
@@ -115 +119 @@
     allsizesy.sort()
     
-    return allsizesx, allsizesy
+    return allsizesx, allsizesy, depression
 
 #########################################################################
@@ -125 +129 @@
 import matplotlib.mlab as mlab
 import mymath as mym
+import myplot as myp
+
+import plfit
+import plplot
+import randht
+import plpva
 
 save = True
 save = False
+pression = False
+pression = True
 
 if save:
-    allsizesx, allsizesy = getsize("psfc.nc")
+    allsizesx, allsizesy, depression = getsize("/home/aymeric/Big_Data/psfc.nc")
     ### sauvegarde texte pour inspection
     mym.writeascii(allsizesx,'allsizex.txt')
     mym.writeascii(allsizesy,'allsizey.txt')
+    mym.writeascii(depression,'alldepression.txt')
     ### sauvegarde binaire pour utilisation python
     myfile = open('allsizex.bin', 'wb')
@@ -141 +154 @@
     pickle.dump(allsizesy, myfile)
     myfile.close()
+    myfile = open('alldepression.bin', 'wb')
+    pickle.dump(depression, myfile)
+    myfile.close()
 
 ### load files
@@ -147 +163 @@
 myfile = open('allsizey.bin', 'r')
 allsizesy = pickle.load(myfile)
-
-### append sizes
-plothist = np.append(allsizesx,allsizesy)
+myfile = open('alldepression.bin', 'r')
+depression = pickle.load(myfile)
+depression = np.array(abs(depression))#*1000.
+
+### sizes
+#plothist = np.append(allsizesx,allsizesy)
 plothist = allsizesx
-#plothist = allsizesy
+if pression: plothist = depression
 plothist.sort()
-
-### make bins
+print 'mean ', np.mean(plothist,dtype=np.float64)  
+print 'std ', np.std(plothist,dtype=np.float64)
+print 'max ', np.max(plothist)
+print 'min ', np.min(plothist)
+print 'len ', len(plothist)
+
+
+### MAKE BINS
 nbins = 100
 zebins = [2.0]
@@ -165 +190 @@
 zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
 nbins = 200
-zebins = [0.2/np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
+
+if pression: zebins = [0.3]
 
 for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
 zebins = np.array(zebins)
-### select reasonable bins for DD
-zebins = zebins [ zebins > 10. ] 
-zebins = zebins [ zebins < 1000. ]
-#print 'bins ',zebins
+#### select reasonable bins for DD
+if not pression:
+    zebins = zebins [ zebins > 15. ] 
+    zebins = zebins [ zebins < 1000. ]
+else:
+    zebins = zebins [ zebins < 10. ]
 print 'corrected bins ',zebins
-print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy) 
-
-#plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.]  ## pour calibrer
 
 #### HISTOGRAM
+plt.figure(1)
 plt.hist(    plothist,\
              log=True,\
              bins=zebins,\
 #             cumulative=-1,\
+             normed=True,\
         )
 plt.xscale('log')
-plt.show()
-
+if pression:    plt.xlabel('Pressure (Pa)')
+else:           plt.xlabel('Diameter (m)')
+plt.ylabel('Population (normalized)')
+if pression: prefix="p"
+else:        prefix=""
+myp.makeplotres(prefix+"histogram",res=200,disp=False)
+plt.close(1)
+
+### COMPARED HISTOGRAMS
+### --- FIT WITH POWER LAW
+if pression: [alpha, xmin, L] = plfit.plfit(plothist,'xmin',0.3)
+else:        [alpha, xmin, L] = plfit.plfit(plothist,'limit',20.)
+print alpha,xmin
+
+#a = plpva.plpva(plothist,0.75,'xmin',0.75)
+#print a
+
+#### DEUXIEME ROUTINE
+####IL FAUT UTILISER LE DISCRET POUR LA TAILLE !!!
+#if pression:   myplfit = plfit.plfit(plothist,verbose=True,xmin=0.75)
+#else:          myplfit = plfit.plfit(plothist,verbose=True,xmin=20.)
+#myplfit.plotppf()
+#plt.show()
+#exit()
+
+
+#plt.figure(1)
+#h = plplot.plplot(plothist,xmin,alpha)
+#myp.makeplotres(prefix+"fit",res=200,disp=False)
+plt.figure(2)
+### --- POWER LAW (factor does not really matter)
+power = (xmin/2.2)*np.array(randht.randht(10000,'powerlaw',alpha))
+#power = (xmin/2.2)*np.array(randht.randht(10000,'cutoff',alpha,10.)) ##marche pas si trop grand
+print 'mean ', np.mean(power,dtype=np.float64)
+### --- EXPONENTIAL LAW
+expo = randht.randht(10000,'exponential',1./(np.mean(power,dtype=np.float64)*1.00))
+print 'mean ', np.mean(expo,dtype=np.float64)
+### --- PLOT
+plt.hist(    [plothist,power,expo],\
+             label=['LES vortices','Power law '+'{:.1f}'.format(alpha),'Exponential law'],\
+             log=True,\
+             bins=zebins,\
+#            cumulative=-1,\
+             normed=True,\
+        )
+plt.legend()
+plt.xscale('log')
+if pression:    plt.xlabel('Pressure (Pa)')
+else:           plt.xlabel('Diameter (m)')
+plt.ylabel('Population (normalized)')
+myp.makeplotres(prefix+"comparison",res=200,disp=False)
+plt.close(2)
+
+########################
+########################
+zebins = [30.,42.,60.,84.,120.,170.,240.,340.]
+plothist = []
+plothist = np.append(plothist,30 *np.ones(306))
+plothist = np.append(plothist,42 *np.ones(58) )
+plothist = np.append(plothist,60 *np.ones(66) )
+plothist = np.append(plothist,84 *np.ones(41) )
+plothist = np.append(plothist,120*np.ones(19) )
+plothist = np.append(plothist,170*np.ones(9)  )
+plothist = np.append(plothist,240*np.ones(2)  )
+plothist = np.append(plothist,340*np.ones(1)  )
+
+#zebins = [50.,71.,100.,141.,200.,282.,400.]
+#plothist = []
+#plothist = np.append(plothist,50. *np.ones(36))
+#plothist = np.append(plothist,71. *np.ones(18) )
+#plothist = np.append(plothist,100. *np.ones(12) )
+#plothist = np.append(plothist,141. *np.ones(6) )
+#plothist = np.append(plothist,200.*np.ones(4) )
+#plothist = np.append(plothist,282.*np.ones(1)  )
+#plothist = np.append(plothist,400.*np.ones(2)  )
+
+plt.figure(3)
+[alpha, xmin, L] = plfit.plfit(plothist,'xmin',30)#50.)
+print alpha,xmin
+#a = plpva.plpva(plothist,30,'xmin',30)
+h = plplot.plplot(plothist,xmin,alpha)
+plt.loglog(h[0], h[1], 'k--',linewidth=2)
+plt.hist(    plothist,\
+             log=True,\
+             bins=zebins,\
+#             cumulative=-1,\
+             normed=True,\
+        )
+plt.xscale('log')
+plt.xlabel('Pressure (micro Pa)')
+plt.ylabel('Population (normalized)')
+myp.makeplotres("data",res=200,disp=False)
+
+#plt.figure(4)
+#[alpha, xmin, L] = plfit.plfit(plothist,'xmin',50.) #,'xmin',30.)
+#print alpha,xmin
+#h = plplot.plplot(plothist,xmin,alpha)
+#myp.makeplotres("datafit",res=200,disp=False)