source: trunk/UTIL/PYTHON/mcd/examples/dustdevil.py @ 1009

#! /usr/bin/env python

from mcd import mcd
from myplot import maplatlon
import numpy as np
import matplotlib.pyplot as mpl

from math import isnan

from ppplot import writeascii


dafile = "/home/aymeric/Work/submitted/coauthorship/2012reiss/reiss_data.txt"

yorgl = np.loadtxt(dafile)
yorgl = np.transpose(yorgl)
ddwind = yorgl[5]
ddwinddir = yorgl[6]
ddwinderror = yorgl[9]

ddheight = yorgl[7]
ddheighterror = yorgl[8]

#ddheight = yorgl[7]*0. + 10.
#ddheighterror = yorgl[8]*0. + 0.5

stddev = np.std(ddwind[np.where(ddwind != -999.)])
print stddev

## missing points are NaN
ind = np.where(ddwind == -999.) ; ddwind[ind] = np.NaN
print "1. wind not available", yorgl[0][ind]
ind = np.where(ddheight == -999.) ; ddwind[ind] = np.NaN
print "1. height not available", yorgl[0][ind]
## no data point with rel error more than x%
#ind = np.where(100. * ddwinderror / ddwinddir > 5.)
#there was an error here!
tol = 50.
ind = np.where(100. * ddwinderror / ddwind > tol) ; ddwind[ind] = np.NaN 
print "2. error on wind too large"
print yorgl[0][ind]
print ddwinderror[ind]
ind = np.where(100. * ddheighterror / ddheight > tol) ; ddwind[ind] = np.NaN 
print "3. error on height too large"
print yorgl[0][ind]
print ddheighterror[ind]
## one record is doubtful according to Reiss
ind = np.where( yorgl[0] == 3 ) ; ddwind[ind] = np.NaN #; ddwinddir[ind] = np.NaN ; ddwinderror[ind] = np.NaN
print "4. doubtful", yorgl[0][ind]

#exit()

## constant_heights
#for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]:
for heights in [10.]:

## multiple_obsheight
#for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]:

 query = mcd()
 mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = [] 

 for i in yorgl[0]:

   if not isnan(ddwind[i-1]):
    query.lat = yorgl[1][i-1]
    query.lon = yorgl[2][i-1]
    query.xdate = yorgl[3][i-1]
    query.loct = yorgl[4][i-1]

    query.xz = ddheight[i-1] #ddheight is really the one that works best

    ## constant_heights
    query.xz = heights

    ## multiple_obsheight
    #query.xz = ddheight[i-1]*heights

    query.update() ; uu = query.zonwind ; vv = query.merwind
    if uu == -999. or vv == -999.: 
       speed = np.NaN ; speedmin = np.NaN ; speedmax = np.NaN
       angle = np.NaN ; anglemin = np.NaN ; anglemax = np.NaN
    else:
       speed = np.sqrt(uu*uu+vv*vv) ; speedmin = speed ; speedmax = speed

       angle = 90. - np.arctan2(vv,uu) * 180. / np.pi
       if angle < 0.: angle = angle + 360.

       #query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
       #query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
       query.xzs = max(query.xz - 0.1*query.xz,10.)
       query.xze = query.xz + 0.1*query.xz

       query.profile()
       tab = np.sqrt(query.zonwindtab*query.zonwindtab + query.merwindtab*query.merwindtab)
       tab2 = 90. - np.arctan2(query.merwindtab,query.zonwindtab) * 180. / np.pi
       tab2[np.where(tab2 < 0.)] = tab2[np.where(tab2 < 0.)] + 360.
       speedmin = min(tab)
       speedmax = max(tab)
       anglemin = min(tab2)-5.
       anglemax = max(tab2)+5.
       if anglemin < 0.: anglemin = angle - (anglemax - angle)
       print int(i), speed, speedmin, speedmax 
       print int(i), angle, anglemin, anglemax
   else:
    #print "removed ",i,ddwinderror[i-1],ddheighterror[i-1]
    speed = -999. ; speedmin = -1000. ; speedmax = -998. ; angle = -999. ; anglemin = -1000. ; anglemax = -998.

   mcdwind.append(speed)
   mcdwindle.append(speed-speedmin)
   mcdwindhe.append(speedmax-speed)
   mcdangle.append(angle)
   mcdanglele.append(angle-anglemin)
   mcdanglehe.append(anglemax-angle)

 mcdwind = np.array(mcdwind)
 mcdwindle = np.array(mcdwindle)
 mcdwindhe = np.array(mcdwindhe)

 mcdangle = np.array(mcdangle)
 mcdanglele = np.array(mcdanglele)
 mcdanglehe = np.array(mcdanglehe)


 mpl.figure(figsize=(12,10))
 mpl.plot([0,30], [stddev,stddev+30], 'g--')
 mpl.plot([0,30], [-stddev,-stddev+30], 'g--')
 mpl.plot([0,30], [0,30], 'g-')

 writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt')


 mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')
# mpl.xlim(xmin=0.,xmax=17.)
# mpl.ylim(ymin=0.,ymax=30.)
 mpl.xlim(xmin=0.,xmax=30.)
 mpl.ylim(ymin=0.,ymax=30.)


 mpl.xlabel('MCD horizontal wind speed (m/s)')
 mpl.ylabel('Dust devil observed drift speed (m/s)')


 for i in yorgl[0]:
    ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord
    ##  qui sont in si on considere l error bar
    if abs(ddwind[i-1] - mcdwind[i-1]) > stddev*1.2:
        mpl.annotate(str(int(i)), xy=(mcdwind[i-1], ddwind[i-1]), xytext=(10,2), 
            textcoords='offset points', ha='center', va='bottom' )#,
            #bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),
            #arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3', 
            #color='red'))


 #mpl.show()
 mpl.savefig("speed"+str(heights)+".eps")
 mpl.savefig("speed"+str(heights)+".png")
 ##mpl.savefig("comp4.eps")
 ##mpl.savefig("comp4.png")

 ddwinddirerror = mcdanglele*0. + 15.

 mpl.figure(figsize=(12,10))
 mpl.errorbar(mcdangle,ddwinddir,yerr=[ddwinddirerror,ddwinddirerror],xerr=[mcdanglele,mcdanglehe],fmt='bo')
 mpl.plot(mcdangle,ddwinddir,'bo')
 mpl.xlim(xmin=0.,xmax=360.)
 mpl.ylim(ymin=0.,ymax=360.)

 mpl.xlabel('MCD horizontal wind orientation (deg)')
 mpl.ylabel('Dust devil observed drift orientation (deg)')


 stddev2 = 50.  ## voir e.g. 8 9 10
 mpl.plot([0,360], [stddev2,stddev2+360], 'g--')
 mpl.plot([0,360], [-stddev2,-stddev2+360], 'g--')
 mpl.plot([0,360], [0,360], 'g-')


 for i in yorgl[0]:
    ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord
    ##  qui sont in si on considere l error bar
    if abs(ddwinddir[i-1] - mcdangle[i-1]) > stddev2*1.2:
        mpl.annotate(str(int(i)), xy=(mcdangle[i-1], ddwinddir[i-1]), xytext=(10,2),
            textcoords='offset points', ha='center', va='bottom')#,
            #bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),
            #arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3',
            #color='red'))

 ##xytext=(18,30)

 #mpl.show()
 mpl.savefig("angle"+str(heights)+".eps")
 mpl.savefig("angle"+str(heights)+".png")
 #mpl.savefig("comp3.eps")
 #mpl.savefig("comp3.png")

#exit()
#
#
#dd = mcd()
#
#
#
#
#
##dd.zkey = 4
##dd.xzs = 600.
##dd.xze = 0.1
##dd.profile()
##dd.plot1d(["t","p"])
##mpl.show()
##exit()
#
### CASE 1
#dd.lat   = -14.6584
#dd.lon   = 175.54
#dd.xdate = 265.278792
#dd.loct  = 14.8487
#dd.xz    = 130.4
#
#dd.update()
#dd.printmcd()
#
#dd.lats  = -5.  
#dd.late  = -25.
#dd.lons  = 160. 
#dd.lone  = 190.
#dd.map2d(["wind"],incwind=True,fixedlt=True)
### il y a un piege pour les cartes locales il faut fixedlt=True
#mpl.savefig("case1.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
#
### CASE 2
#dd.lat   = -68.6125
#dd.lon   = 11.4303
#dd.xdate = 286.507293
#dd.loct  = 15.129030
#dd.xz    = 872.8
##dd.xz    = 10000.0
#
#dd.update()
#dd.printmcd()
#
#dd.lats  = -75.
#dd.late  = -60.
#dd.lons  = 0.
#dd.lone  = 20.
#dd.map2d(["wind"],incwind=True,fixedlt=True)
#mpl.savefig("case2.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
#
###3    25,1306 314,4842        60,8946 14,9971 30,89575722     232     250,217961736111
#
### CASE 3
#dd.lat   = 25.1306
#dd.lon   = 314.4842
#dd.xdate = 60.8946
#dd.loct  = 14.9971
#dd.xz    = 250.218
#
#dd.update()
#dd.printmcd()
#
###4    35,564  199,486 60,5883 14,9949 12,2972688101353        125     895,247012611329
#
### CASE 4
#dd.lat   = 35.564
#dd.lon   = 199.486
#dd.xdate = 60.5883
#dd.loct  = 14.9949
#dd.xz    = 895.247
#
#dd.update()
#dd.printmcd()
#
###5    68,346  234,396 142,828 15,1777 13,4079899322222        83      1128,06581216829
###6    68,316  234,462 142,828 15,1819 16,1939071396022        85      582,624074786015
###7    68,302  234,144 142,828 15,1607 17,3354121022885        112     262,299040101764
#
### CASE 5
#dd.lat   = 68.32
#dd.lon   = 234.3
#dd.xdate = 142.828
#dd.loct  = 15.2
#
#dd.xz    = 1128.066
#dd.update()
#dd.printmcd()
#
#dd.xz    = 582.624
#dd.update()
#dd.printmcd()
#
#dd.xz    = 262.299
#dd.update()
#dd.printmcd()
#
#
###19   27,055  129,614 13,5818 14,444  42,6488205391137        302     1395,96076100437
#
### CASE 19
#dd.lat   = 27.055
#dd.lon   = 129.614
#dd.xdate = 13.5818
#dd.loct  = 14.444
#
#dd.xz    = 1395.961
#dd.update()
#dd.printmcd()
#
#dd.lats  = 20.
#dd.late  = 30.
#dd.lons  = 125.
#dd.lone  = 135.
#dd.map2d(["wind"],incwind=True,fixedlt=True)
#mpl.savefig("case19.png",dpi=110,bbox_inches='tight',pad_inches=0.4)
#
#