Changeset 1007 for trunk/UTIL


Ignore:
Timestamp:
Jul 19, 2013, 10:35:53 AM (11 years ago)
Author:
aslmd
Message:

UTIL PYTHON. various updates on example + mcd + plot scripts. nothing major.

Location:
trunk/UTIL/PYTHON
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/mcd/examples/dustdevil.py

    r943 r1007  
    77
    88from math import isnan
     9
     10from ppplot import writeascii
     11
    912
    1013dafile = "/home/aymeric/Work/submitted/coauthorship/2012reiss/reiss_data.txt"
    … …  
    4851#exit()
    4952
    50 query = mcd()
    51 mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = []
    52 
    53 for i in yorgl[0]:
     53## constant_heights
     54#for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]:
     55for heights in [10.]:
     56
     57## multiple_obsheight
     58#for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]:
     59
     60 query = mcd()
     61 mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = []
     62
     63 for i in yorgl[0]:
    5464
    5565   if not isnan(ddwind[i-1]):
    … …  
    5868    query.xdate = yorgl[3][i-1]
    5969    query.loct = yorgl[4][i-1]
    60     query.xz = ddheight[i-1]
    61     #query.xz = 1000.  #ddheight is really the one that works best
     70
     71    query.xz = ddheight[i-1] #ddheight is really the one that works best
     72
     73    ## constant_heights
     74    query.xz = heights
     75
     76    ## multiple_obsheight
     77    #query.xz = ddheight[i-1]*heights
     78
    6279    query.update() ; uu = query.zonwind ; vv = query.merwind
    6380    if uu == -999. or vv == -999.:
    … …  
    7087       if angle < 0.: angle = angle + 360.
    7188
    72        query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
    73        query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
     89       #query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
     90       #query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
     91       query.xzs = max(query.xz - 0.1*query.xz,10.)
     92       query.xze = query.xz + 0.1*query.xz
     93
    7494       query.profile()
    7595       tab = np.sqrt(query.zonwindtab*query.zonwindtab + query.merwindtab*query.merwindtab)
    … …  
    87107    speed = -999. ; speedmin = -1000. ; speedmax = -998. ; angle = -999. ; anglemin = -1000. ; anglemax = -998.
    88108
    89    mcdwind = np.append(mcdwind,speed)
    90    mcdwindle = np.append(mcdwindle,speed-speedmin)
    91    mcdwindhe = np.append(mcdwindhe,speedmax-speed)
    92    mcdangle = np.append(mcdangle,angle)
    93    mcdanglele = np.append(mcdanglele,angle-anglemin)
    94    mcdanglehe = np.append(mcdanglehe,anglemax-angle)
    95 
    96 mpl.figure(figsize=(12,10))
    97 mpl.plot([0,30], [stddev,stddev+30], 'g--')
    98 mpl.plot([0,30], [-stddev,-stddev+30], 'g--')
    99 mpl.plot([0,30], [0,30], 'g-')
    100 
    101 mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')
    102 mpl.xlim(xmin=0.,xmax=17.)
    103 mpl.ylim(ymin=0.,ymax=30.)
    104 
    105 mpl.xlabel('MCD horizontal wind speed (m/s)')
    106 mpl.ylabel('Dust devil observed drift speed (m/s)')
    107 
    108 
    109 for i in yorgl[0]:
     109   mcdwind.append(speed)
     110   mcdwindle.append(speed-speedmin)
     111   mcdwindhe.append(speedmax-speed)
     112   mcdangle.append(angle)
     113   mcdanglele.append(angle-anglemin)
     114   mcdanglehe.append(anglemax-angle)
     115
     116 mcdwind = np.array(mcdwind)
     117 mcdwindle = np.array(mcdwindle)
     118 mcdwindhe = np.array(mcdwindhe)
     119
     120 mcdangle = np.array(mcdangle)
     121 mcdanglele = np.array(mcdanglele)
     122 mcdanglehe = np.array(mcdanglehe)
     123
     124
     125 mpl.figure(figsize=(12,10))
     126 mpl.plot([0,30], [stddev,stddev+30], 'g--')
     127 mpl.plot([0,30], [-stddev,-stddev+30], 'g--')
     128 mpl.plot([0,30], [0,30], 'g-')
     129
     130 writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt')
     131
     132
     133 mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')
     134# mpl.xlim(xmin=0.,xmax=17.)
     135# mpl.ylim(ymin=0.,ymax=30.)
     136 mpl.xlim(xmin=0.,xmax=30.)
     137 mpl.ylim(ymin=0.,ymax=30.)
     138
     139
     140 mpl.xlabel('MCD horizontal wind speed (m/s)')
     141 mpl.ylabel('Dust devil observed drift speed (m/s)')
     142
     143
     144 for i in yorgl[0]:
    110145    ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord
    111146    ##  qui sont in si on considere l error bar
    … …  
    118153
    119154
    120 ##mpl.show()
    121 mpl.savefig("comp.eps")
    122 mpl.savefig("comp.png")
    123 #mpl.savefig("comp4.eps")
    124 #mpl.savefig("comp4.png")
    125 
    126 ddwinddirerror = mcdanglele*0. + 15.
    127 
    128 mpl.figure(figsize=(12,10))
    129 mpl.errorbar(mcdangle,ddwinddir,yerr=[ddwinddirerror,ddwinddirerror],xerr=[mcdanglele,mcdanglehe],fmt='bo')
    130 mpl.plot(mcdangle,ddwinddir,'bo')
    131 mpl.xlim(xmin=0.,xmax=360.)
    132 mpl.ylim(ymin=0.,ymax=360.)
    133 
    134 mpl.xlabel('MCD horizontal wind orientation (deg)')
    135 mpl.ylabel('Dust devil observed drift orientation (deg)')
    136 
    137 
    138 stddev = 50.  ## voir e.g. 8 9 10
    139 mpl.plot([0,360], [stddev,stddev+360], 'g--')
    140 mpl.plot([0,360], [-stddev,-stddev+360], 'g--')
    141 mpl.plot([0,360], [0,360], 'g-')
    142 
    143 
    144 for i in yorgl[0]:
     155 #mpl.show()
     156 mpl.savefig("speed"+str(heights)+".eps")
     157 mpl.savefig("speed"+str(heights)+".png")
     158 ##mpl.savefig("comp4.eps")
     159 ##mpl.savefig("comp4.png")
     160
     161 ddwinddirerror = mcdanglele*0. + 15.
     162
     163 mpl.figure(figsize=(12,10))
     164 mpl.errorbar(mcdangle,ddwinddir,yerr=[ddwinddirerror,ddwinddirerror],xerr=[mcdanglele,mcdanglehe],fmt='bo')
     165 mpl.plot(mcdangle,ddwinddir,'bo')
     166 mpl.xlim(xmin=0.,xmax=360.)
     167 mpl.ylim(ymin=0.,ymax=360.)
     168
     169 mpl.xlabel('MCD horizontal wind orientation (deg)')
     170 mpl.ylabel('Dust devil observed drift orientation (deg)')
     171
     172
     173 stddev2 = 50.  ## voir e.g. 8 9 10
     174 mpl.plot([0,360], [stddev2,stddev2+360], 'g--')
     175 mpl.plot([0,360], [-stddev2,-stddev2+360], 'g--')
     176 mpl.plot([0,360], [0,360], 'g-')
     177
     178
     179 for i in yorgl[0]:
    145180    ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord
    146181    ##  qui sont in si on considere l error bar
    147     if abs(ddwinddir[i-1] - mcdangle[i-1]) > stddev*1.2:
     182    if abs(ddwinddir[i-1] - mcdangle[i-1]) > stddev2*1.2:
    148183        mpl.annotate(str(int(i)), xy=(mcdangle[i-1], ddwinddir[i-1]), xytext=(10,2),
    149184            textcoords='offset points', ha='center', va='bottom')#,
    … …  
    152187            #color='red'))
    153188
    154 ##xytext=(18,30)
    155 
    156 
    157 mpl.savefig("comp2.eps")
    158 mpl.savefig("comp2.png")
    159 #mpl.savefig("comp3.eps")
    160 #mpl.savefig("comp3.png")
    161 
    162 
    163 exit()
    164 
    165 
    166 dd = mcd()
    167 
    168 
    169 
    170 
    171 
    172 #dd.zkey = 4
    173 #dd.xzs = 600.
    174 #dd.xze = 0.1
    175 #dd.profile()
    176 #dd.plot1d(["t","p"])
    177 #mpl.show()
     189 ##xytext=(18,30)
     190
     191 #mpl.show()
     192 mpl.savefig("angle"+str(heights)+".eps")
     193 mpl.savefig("angle"+str(heights)+".png")
     194 #mpl.savefig("comp3.eps")
     195 #mpl.savefig("comp3.png")
     196
    178197#exit()
    179 
    180 ## CASE 1
    181 dd.lat   = -14.6584
    182 dd.lon   = 175.54
    183 dd.xdate = 265.278792
    184 dd.loct  = 14.8487
    185 dd.xz    = 130.4
    186 
    187 dd.update()
    188 dd.printmcd()
    189 
    190 dd.lats  = -5. 
    191 dd.late  = -25.
    192 dd.lons  = 160.
    193 dd.lone  = 190.
    194 dd.map2d(["wind"],incwind=True,fixedlt=True)
    195 ## il y a un piege pour les cartes locales il faut fixedlt=True
    196 mpl.savefig("case1.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
    197 
    198 ## CASE 2
    199 dd.lat   = -68.6125
    200 dd.lon   = 11.4303
    201 dd.xdate = 286.507293
    202 dd.loct  = 15.129030
    203 dd.xz    = 872.8
    204 #dd.xz    = 10000.0
    205 
    206 dd.update()
    207 dd.printmcd()
    208 
    209 dd.lats  = -75.
    210 dd.late  = -60.
    211 dd.lons  = 0.
    212 dd.lone  = 20.
    213 dd.map2d(["wind"],incwind=True,fixedlt=True)
    214 mpl.savefig("case2.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
    215 
    216 ##3     25,1306 314,4842        60,8946 14,9971 30,89575722     232     250,217961736111
    217 
    218 ## CASE 3
    219 dd.lat   = 25.1306
    220 dd.lon   = 314.4842
    221 dd.xdate = 60.8946
    222 dd.loct  = 14.9971
    223 dd.xz    = 250.218
    224 
    225 dd.update()
    226 dd.printmcd()
    227 
    228 ##4     35,564  199,486 60,5883 14,9949 12,2972688101353        125     895,247012611329
    229 
    230 ## CASE 4
    231 dd.lat   = 35.564
    232 dd.lon   = 199.486
    233 dd.xdate = 60.5883
    234 dd.loct  = 14.9949
    235 dd.xz    = 895.247
    236 
    237 dd.update()
    238 dd.printmcd()
    239 
    240 ##5     68,346  234,396 142,828 15,1777 13,4079899322222        83      1128,06581216829
    241 ##6     68,316  234,462 142,828 15,1819 16,1939071396022        85      582,624074786015
    242 ##7     68,302  234,144 142,828 15,1607 17,3354121022885        112     262,299040101764
    243 
    244 ## CASE 5
    245 dd.lat   = 68.32
    246 dd.lon   = 234.3
    247 dd.xdate = 142.828
    248 dd.loct  = 15.2
    249 
    250 dd.xz    = 1128.066
    251 dd.update()
    252 dd.printmcd()
    253 
    254 dd.xz    = 582.624
    255 dd.update()
    256 dd.printmcd()
    257 
    258 dd.xz    = 262.299
    259 dd.update()
    260 dd.printmcd()
    261 
    262 
    263 ##19    27,055  129,614 13,5818 14,444  42,6488205391137        302     1395,96076100437
    264 
    265 ## CASE 19
    266 dd.lat   = 27.055
    267 dd.lon   = 129.614
    268 dd.xdate = 13.5818
    269 dd.loct  = 14.444
    270 
    271 dd.xz    = 1395.961
    272 dd.update()
    273 dd.printmcd()
    274 
    275 dd.lats  = 20.
    276 dd.late  = 30.
    277 dd.lons  = 125.
    278 dd.lone  = 135.
    279 dd.map2d(["wind"],incwind=True,fixedlt=True)
    280 mpl.savefig("case19.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
    281 
    282 
     198#
     199#
     200#dd = mcd()
     201#
     202#
     203#
     204#
     205#
     206##dd.zkey = 4
     207##dd.xzs = 600.
     208##dd.xze = 0.1
     209##dd.profile()
     210##dd.plot1d(["t","p"])
     211##mpl.show()
     212##exit()
     213#
     214### CASE 1
     215#dd.lat   = -14.6584
     216#dd.lon   = 175.54
     217#dd.xdate = 265.278792
     218#dd.loct  = 14.8487
     219#dd.xz    = 130.4
     220#
     221#dd.update()
     222#dd.printmcd()
     223#
     224#dd.lats  = -5. 
     225#dd.late  = -25.
     226#dd.lons  = 160.
     227#dd.lone  = 190.
     228#dd.map2d(["wind"],incwind=True,fixedlt=True)
     229### il y a un piege pour les cartes locales il faut fixedlt=True
     230#mpl.savefig("case1.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
     231#
     232### CASE 2
     233#dd.lat   = -68.6125
     234#dd.lon   = 11.4303
     235#dd.xdate = 286.507293
     236#dd.loct  = 15.129030
     237#dd.xz    = 872.8
     238##dd.xz    = 10000.0
     239#
     240#dd.update()
     241#dd.printmcd()
     242#
     243#dd.lats  = -75.
     244#dd.late  = -60.
     245#dd.lons  = 0.
     246#dd.lone  = 20.
     247#dd.map2d(["wind"],incwind=True,fixedlt=True)
     248#mpl.savefig("case2.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
     249#
     250###3    25,1306 314,4842        60,8946 14,9971 30,89575722     232     250,217961736111
     251#
     252### CASE 3
     253#dd.lat   = 25.1306
     254#dd.lon   = 314.4842
     255#dd.xdate = 60.8946
     256#dd.loct  = 14.9971
     257#dd.xz    = 250.218
     258#
     259#dd.update()
     260#dd.printmcd()
     261#
     262###4    35,564  199,486 60,5883 14,9949 12,2972688101353        125     895,247012611329
     263#
     264### CASE 4
     265#dd.lat   = 35.564
     266#dd.lon   = 199.486
     267#dd.xdate = 60.5883
     268#dd.loct  = 14.9949
     269#dd.xz    = 895.247
     270#
     271#dd.update()
     272#dd.printmcd()
     273#
     274###5    68,346  234,396 142,828 15,1777 13,4079899322222        83      1128,06581216829
     275###6    68,316  234,462 142,828 15,1819 16,1939071396022        85      582,624074786015
     276###7    68,302  234,144 142,828 15,1607 17,3354121022885        112     262,299040101764
     277#
     278### CASE 5
     279#dd.lat   = 68.32
     280#dd.lon   = 234.3
     281#dd.xdate = 142.828
     282#dd.loct  = 15.2
     283#
     284#dd.xz    = 1128.066
     285#dd.update()
     286#dd.printmcd()
     287#
     288#dd.xz    = 582.624
     289#dd.update()
     290#dd.printmcd()
     291#
     292#dd.xz    = 262.299
     293#dd.update()
     294#dd.printmcd()
     295#
     296#
     297###19   27,055  129,614 13,5818 14,444  42,6488205391137        302     1395,96076100437
     298#
     299### CASE 19
     300#dd.lat   = 27.055
     301#dd.lon   = 129.614
     302#dd.xdate = 13.5818
     303#dd.loct  = 14.444
     304#
     305#dd.xz    = 1395.961
     306#dd.update()
     307#dd.printmcd()
     308#
     309#dd.lats  = 20.
     310#dd.late  = 30.
     311#dd.lons  = 125.
     312#dd.lone  = 135.
     313#dd.map2d(["wind"],incwind=True,fixedlt=True)
     314#mpl.savefig("case19.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
     315#
     316#

  • trunk/UTIL/PYTHON/mcd/examples/gale.py

    r943 r1007  
    3535mpl.savefig("temp.png",dpi=85,bbox_inches='tight',pad_inches=0.25)
    3636
     37exit()
     38
     39
     40
    3741
    3842#gale.xdate = 270.

  • trunk/UTIL/PYTHON/mcd/inimeso/inimeso.py

    r943 r1007  
    3737sounding.close() ; additional.close()
    3838query.plot1d(["p","t","u","v"]) ; mpl.show()
     39
     40### Add information about dust opacity
     41dod = open("dustopacity.def", "w")
     42dod.write( "%4.2f"  %(query.extvartab[iz,36]) )
     43dod.close()

  • trunk/UTIL/PYTHON/planetoplot_v2/ppcompute.py

    r1002 r1007  
    9595           else:
    9696              return np.array(field).sum(axis=axis)
     97
     98## compute mean over bins
     99## author AS
     100def meanbin(y,x,bins):
     101    bins.sort() # sorting is needed for binning
     102    meanvalues = []
     103    for iii in range(len(bins)):
     104        ## GET VALUES OVER RELEVANT INTERVALS
     105        if iii == 0:
     106          ind = x<bins[iii+1] ; ys = y[ind]
     107        elif iii == len(bins)-1:
     108          ind = x>=bins[iii-1] ; ys = y[ind]
     109        else:
     110          ind = x>=bins[iii-1] ; interm = x[ind] ; intermf = y[ind]
     111          ind = interm<bins[iii+1] ; xs = interm[ind] ; ys = intermf[ind]
     112        ## COMPUTE MEAN and TREAT NAN CASE
     113        meanvalues.append(mean(ys))
     114    ## RETURN A NUMPY ARRAY
     115    meanvalues = np.array(meanvalues)
     116    return meanvalues
    97117
    98118################

  • trunk/UTIL/PYTHON/planetoplot_v2/ppplot.py

    r1004 r1007  
    180180# useful functions #
    181181####################
     182
     183# continuity with matplotlib
     184def close():
     185    mpl.close()
     186
     187# a function for predefined figure sizes
     188def figuref(x=16,y=9):
     189    mpl.figure(figsize=(x,y))
    182190
    183191# a function to change color lines according to a color map (idea by A. Pottier)
    … …  
    649657                objC = mpl.contour(x, y, what_I_contour, \
    650658                            zelevelsc, colors = ccol, linewidths = cline)
    651                 #mpl.clabel(objC, inline=1, fontsize=10)
     659                #mpl.clabel(objC, inline=1, fontsize="small",\
     660                #             inline_spacing=1,fmt="%.0f")
    652661            mpl.contourf(x, y, \
    653662                         self.field, \
    654663                         zelevels, cmap=palette)
     664            #mpl.pcolor(x,y,\
     665            #             self.field, \
     666            #             cmap=palette)
    655667            # make log axes and/or invert ordinate
    656668            ax = mpl.gca()
Note: See TracChangeset for help on using the changeset viewer.