| 1 | #! /usr/bin/env python |
|---|
| 2 | |
|---|
| 3 | from mcd import mcd |
|---|
| 4 | from myplot import maplatlon |
|---|
| 5 | import numpy as np |
|---|
| 6 | import matplotlib.pyplot as mpl |
|---|
| 7 | |
|---|
| 8 | from math import isnan |
|---|
| 9 | |
|---|
| 10 | from ppplot import writeascii |
|---|
| 11 | |
|---|
| 12 | |
|---|
| 13 | from scipy import stats |
|---|
| 14 | |
|---|
| 15 | |
|---|
| 16 | dafile = "/home/aymeric/Work/submitted/coauthorship/2012reiss/reiss_data.txt" |
|---|
| 17 | |
|---|
| 18 | yorgl = np.loadtxt(dafile) |
|---|
| 19 | yorgl = np.transpose(yorgl) |
|---|
| 20 | ddwind = yorgl[5] |
|---|
| 21 | ddwinddir = yorgl[6] |
|---|
| 22 | ddwinderror = yorgl[9] |
|---|
| 23 | |
|---|
| 24 | ddheight = yorgl[7] |
|---|
| 25 | ddheighterror = yorgl[8] |
|---|
| 26 | |
|---|
| 27 | #ddheight = yorgl[7]*0. + 10. |
|---|
| 28 | #ddheighterror = yorgl[8]*0. + 0.5 |
|---|
| 29 | |
|---|
| 30 | stddev = np.std(ddwind[np.where(ddwind != -999.)]) |
|---|
| 31 | print stddev |
|---|
| 32 | |
|---|
| 33 | ## missing points are NaN |
|---|
| 34 | ind = np.where(ddwind == -999.) ; ddwind[ind] = np.NaN |
|---|
| 35 | print "1. wind not available", yorgl[0][ind] |
|---|
| 36 | ind = np.where(ddheight == -999.) ; ddwind[ind] = np.NaN |
|---|
| 37 | print "1. height not available", yorgl[0][ind] |
|---|
| 38 | ## no data point with rel error more than x% |
|---|
| 39 | #ind = np.where(100. * ddwinderror / ddwinddir > 5.) |
|---|
| 40 | #there was an error here! |
|---|
| 41 | tol = 50. |
|---|
| 42 | ind = np.where(100. * ddwinderror / ddwind > tol) ; ddwind[ind] = np.NaN |
|---|
| 43 | print "2. error on wind too large" |
|---|
| 44 | print yorgl[0][ind] |
|---|
| 45 | print ddwinderror[ind] |
|---|
| 46 | ind = np.where(100. * ddheighterror / ddheight > tol) ; ddwind[ind] = np.NaN |
|---|
| 47 | print "3. error on height too large" |
|---|
| 48 | print yorgl[0][ind] |
|---|
| 49 | print ddheighterror[ind] |
|---|
| 50 | ## one record is doubtful according to Reiss |
|---|
| 51 | ind = np.where( yorgl[0] == 3 ) ; ddwind[ind] = np.NaN #; ddwinddir[ind] = np.NaN ; ddwinderror[ind] = np.NaN |
|---|
| 52 | print "4. doubtful", yorgl[0][ind] |
|---|
| 53 | |
|---|
| 54 | #exit() |
|---|
| 55 | |
|---|
| 56 | ## constant_heights |
|---|
| 57 | #for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]: |
|---|
| 58 | #for heights in [10.]: |
|---|
| 59 | |
|---|
| 60 | ## multiple_obsheight |
|---|
| 61 | #for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]: |
|---|
| 62 | for heights in [1.0]: |
|---|
| 63 | |
|---|
| 64 | query = mcd() |
|---|
| 65 | mcdwind = [] ; mcdwindle = [] ; mcdwindhe = [] ; mcdangle = [] ; mcdanglele = [] ; mcdanglehe = [] |
|---|
| 66 | |
|---|
| 67 | for i in yorgl[0]: |
|---|
| 68 | |
|---|
| 69 | if not isnan(ddwind[i-1]): |
|---|
| 70 | query.lat = yorgl[1][i-1] |
|---|
| 71 | query.lon = yorgl[2][i-1] |
|---|
| 72 | query.xdate = yorgl[3][i-1] |
|---|
| 73 | query.loct = yorgl[4][i-1] |
|---|
| 74 | |
|---|
| 75 | query.xz = ddheight[i-1] #ddheight is really the one that works best |
|---|
| 76 | |
|---|
| 77 | ## constant_heights |
|---|
| 78 | #query.xz = heights |
|---|
| 79 | ## multiple_obsheight |
|---|
| 80 | query.xz = ddheight[i-1]*heights |
|---|
| 81 | |
|---|
| 82 | query.update() ; uu = query.zonwind ; vv = query.merwind |
|---|
| 83 | if uu == -999. or vv == -999.: |
|---|
| 84 | speed = np.NaN ; speedmin = np.NaN ; speedmax = np.NaN |
|---|
| 85 | angle = np.NaN ; anglemin = np.NaN ; anglemax = np.NaN |
|---|
| 86 | else: |
|---|
| 87 | speed = np.sqrt(uu*uu+vv*vv) ; speedmin = speed ; speedmax = speed |
|---|
| 88 | |
|---|
| 89 | angle = 90. - np.arctan2(vv,uu) * 180. / np.pi |
|---|
| 90 | if angle < 0.: angle = angle + 360. |
|---|
| 91 | |
|---|
| 92 | ## multiple_obsheight |
|---|
| 93 | query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.) |
|---|
| 94 | query.xze = ddheight[i-1] + 2.*ddheighterror[i-1] |
|---|
| 95 | ## constant_heights |
|---|
| 96 | #query.xzs = max(query.xz - 0.1*query.xz,10.) |
|---|
| 97 | #query.xze = query.xz + 0.1*query.xz |
|---|
| 98 | |
|---|
| 99 | query.profile() |
|---|
| 100 | tab = np.sqrt(query.zonwindtab*query.zonwindtab + query.merwindtab*query.merwindtab) |
|---|
| 101 | tab2 = 90. - np.arctan2(query.merwindtab,query.zonwindtab) * 180. / np.pi |
|---|
| 102 | tab2[np.where(tab2 < 0.)] = tab2[np.where(tab2 < 0.)] + 360. |
|---|
| 103 | speedmin = min(tab) |
|---|
| 104 | speedmax = max(tab) |
|---|
| 105 | anglemin = min(tab2)-5. |
|---|
| 106 | anglemax = max(tab2)+5. |
|---|
| 107 | if anglemin < 0.: anglemin = angle - (anglemax - angle) |
|---|
| 108 | print int(i), speed, speedmin, speedmax |
|---|
| 109 | print int(i), angle, anglemin, anglemax |
|---|
| 110 | else: |
|---|
| 111 | #print "removed ",i,ddwinderror[i-1],ddheighterror[i-1] |
|---|
| 112 | speed = -999. ; speedmin = -1000. ; speedmax = -998. ; angle = -999. ; anglemin = -1000. ; anglemax = -998. |
|---|
| 113 | |
|---|
| 114 | mcdwind.append(speed) |
|---|
| 115 | mcdwindle.append(speed-speedmin) |
|---|
| 116 | mcdwindhe.append(speedmax-speed) |
|---|
| 117 | mcdangle.append(angle) |
|---|
| 118 | mcdanglele.append(angle-anglemin) |
|---|
| 119 | mcdanglehe.append(anglemax-angle) |
|---|
| 120 | |
|---|
| 121 | mcdwind = np.array(mcdwind) |
|---|
| 122 | mcdwindle = np.array(mcdwindle) |
|---|
| 123 | mcdwindhe = np.array(mcdwindhe) |
|---|
| 124 | |
|---|
| 125 | mcdangle = np.array(mcdangle) |
|---|
| 126 | mcdanglele = np.array(mcdanglele) |
|---|
| 127 | mcdanglehe = np.array(mcdanglehe) |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | mpl.figure(figsize=(12,10)) |
|---|
| 131 | mpl.plot([0,30], [stddev,stddev+30], 'g--') |
|---|
| 132 | mpl.plot([0,30], [-stddev,-stddev+30], 'g--') |
|---|
| 133 | mpl.plot([0,30], [0,30], 'g-') |
|---|
| 134 | |
|---|
| 135 | writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt') |
|---|
| 136 | |
|---|
| 137 | #ind = np.isnan(ddwind) |
|---|
| 138 | #y = ddwind |
|---|
| 139 | #y[ind] = -9999. |
|---|
| 140 | #xi = mcdwind[y > -9999.] |
|---|
| 141 | #y = ddwind[y > -9999.] |
|---|
| 142 | #slope, intercept, r_value, p_value, std_err = stats.linregress(xi,y) |
|---|
| 143 | #linear = slope*xi+intercept |
|---|
| 144 | #mpl.plot(xi,linear) |
|---|
| 145 | #ecart = y - linear |
|---|
| 146 | #print np.std(ecart) |
|---|
| 147 | |
|---|
| 148 | mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo') |
|---|
| 149 | # mpl.xlim(xmin=0.,xmax=17.) |
|---|
| 150 | # mpl.ylim(ymin=0.,ymax=30.) |
|---|
| 151 | mpl.xlim(xmin=0.,xmax=30.) |
|---|
| 152 | mpl.ylim(ymin=0.,ymax=30.) |
|---|
| 153 | |
|---|
| 154 | |
|---|
| 155 | mpl.xlabel('MCD horizontal wind speed (m/s)') |
|---|
| 156 | mpl.ylabel('Dust devil observed drift speed (m/s)') |
|---|
| 157 | |
|---|
| 158 | |
|---|
| 159 | for i in yorgl[0]: |
|---|
| 160 | ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord |
|---|
| 161 | ## qui sont in si on considere l error bar |
|---|
| 162 | if abs(ddwind[i-1] - mcdwind[i-1]) > stddev*1.2: |
|---|
| 163 | mpl.annotate(str(int(i)), xy=(mcdwind[i-1], ddwind[i-1]), xytext=(10,2), |
|---|
| 164 | textcoords='offset points', ha='center', va='bottom' )#, |
|---|
| 165 | #bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3), |
|---|
| 166 | #arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3', |
|---|
| 167 | #color='red')) |
|---|
| 168 | |
|---|
| 169 | |
|---|
| 170 | #mpl.show() |
|---|
| 171 | mpl.savefig("speed"+str(heights)+".eps") |
|---|
| 172 | mpl.savefig("speed"+str(heights)+".png") |
|---|
| 173 | ##mpl.savefig("comp4.eps") |
|---|
| 174 | ##mpl.savefig("comp4.png") |
|---|
| 175 | |
|---|
| 176 | ddwinddirerror = mcdanglele*0. + 15. |
|---|
| 177 | |
|---|
| 178 | mpl.figure(figsize=(12,10)) |
|---|
| 179 | mpl.errorbar(mcdangle,ddwinddir,yerr=[ddwinddirerror,ddwinddirerror],xerr=[mcdanglele,mcdanglehe],fmt='bo') |
|---|
| 180 | mpl.plot(mcdangle,ddwinddir,'bo') |
|---|
| 181 | mpl.xlim(xmin=0.,xmax=360.) |
|---|
| 182 | mpl.ylim(ymin=0.,ymax=360.) |
|---|
| 183 | |
|---|
| 184 | mpl.xlabel('MCD horizontal wind orientation (deg)') |
|---|
| 185 | mpl.ylabel('Dust devil observed drift orientation (deg)') |
|---|
| 186 | |
|---|
| 187 | |
|---|
| 188 | stddev2 = 50. ## voir e.g. 8 9 10 |
|---|
| 189 | mpl.plot([0,360], [stddev2,stddev2+360], 'g--') |
|---|
| 190 | mpl.plot([0,360], [-stddev2,-stddev2+360], 'g--') |
|---|
| 191 | mpl.plot([0,360], [0,360], 'g-') |
|---|
| 192 | |
|---|
| 193 | |
|---|
| 194 | for i in yorgl[0]: |
|---|
| 195 | ## on multiplie stddev par 1.2 de facon a enlever les points tres au bord |
|---|
| 196 | ## qui sont in si on considere l error bar |
|---|
| 197 | if abs(ddwinddir[i-1] - mcdangle[i-1]) > stddev2*1.2: |
|---|
| 198 | mpl.annotate(str(int(i)), xy=(mcdangle[i-1], ddwinddir[i-1]), xytext=(10,2), |
|---|
| 199 | textcoords='offset points', ha='center', va='bottom')#, |
|---|
| 200 | #bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3), |
|---|
| 201 | #arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3', |
|---|
| 202 | #color='red')) |
|---|
| 203 | |
|---|
| 204 | ##xytext=(18,30) |
|---|
| 205 | |
|---|
| 206 | #mpl.show() |
|---|
| 207 | mpl.savefig("angle"+str(heights)+".eps") |
|---|
| 208 | mpl.savefig("angle"+str(heights)+".png") |
|---|
| 209 | #mpl.savefig("comp3.eps") |
|---|
| 210 | #mpl.savefig("comp3.png") |
|---|
| 211 | |
|---|
| 212 | #exit() |
|---|
| 213 | # |
|---|
| 214 | # |
|---|
| 215 | #dd = mcd() |
|---|
| 216 | # |
|---|
| 217 | # |
|---|
| 218 | # |
|---|
| 219 | # |
|---|
| 220 | # |
|---|
| 221 | ##dd.zkey = 4 |
|---|
| 222 | ##dd.xzs = 600. |
|---|
| 223 | ##dd.xze = 0.1 |
|---|
| 224 | ##dd.profile() |
|---|
| 225 | ##dd.plot1d(["t","p"]) |
|---|
| 226 | ##mpl.show() |
|---|
| 227 | ##exit() |
|---|
| 228 | # |
|---|
| 229 | ### CASE 1 |
|---|
| 230 | #dd.lat = -14.6584 |
|---|
| 231 | #dd.lon = 175.54 |
|---|
| 232 | #dd.xdate = 265.278792 |
|---|
| 233 | #dd.loct = 14.8487 |
|---|
| 234 | #dd.xz = 130.4 |
|---|
| 235 | # |
|---|
| 236 | #dd.update() |
|---|
| 237 | #dd.printmcd() |
|---|
| 238 | # |
|---|
| 239 | #dd.lats = -5. |
|---|
| 240 | #dd.late = -25. |
|---|
| 241 | #dd.lons = 160. |
|---|
| 242 | #dd.lone = 190. |
|---|
| 243 | #dd.map2d(["wind"],incwind=True,fixedlt=True) |
|---|
| 244 | ### il y a un piege pour les cartes locales il faut fixedlt=True |
|---|
| 245 | #mpl.savefig("case1.png",dpi=110,bbox_inches='tight',pad_inches=0.4) |
|---|
| 246 | # |
|---|
| 247 | ### CASE 2 |
|---|
| 248 | #dd.lat = -68.6125 |
|---|
| 249 | #dd.lon = 11.4303 |
|---|
| 250 | #dd.xdate = 286.507293 |
|---|
| 251 | #dd.loct = 15.129030 |
|---|
| 252 | #dd.xz = 872.8 |
|---|
| 253 | ##dd.xz = 10000.0 |
|---|
| 254 | # |
|---|
| 255 | #dd.update() |
|---|
| 256 | #dd.printmcd() |
|---|
| 257 | # |
|---|
| 258 | #dd.lats = -75. |
|---|
| 259 | #dd.late = -60. |
|---|
| 260 | #dd.lons = 0. |
|---|
| 261 | #dd.lone = 20. |
|---|
| 262 | #dd.map2d(["wind"],incwind=True,fixedlt=True) |
|---|
| 263 | #mpl.savefig("case2.png",dpi=110,bbox_inches='tight',pad_inches=0.4) |
|---|
| 264 | # |
|---|
| 265 | ###3 25,1306 314,4842 60,8946 14,9971 30,89575722 232 250,217961736111 |
|---|
| 266 | # |
|---|
| 267 | ### CASE 3 |
|---|
| 268 | #dd.lat = 25.1306 |
|---|
| 269 | #dd.lon = 314.4842 |
|---|
| 270 | #dd.xdate = 60.8946 |
|---|
| 271 | #dd.loct = 14.9971 |
|---|
| 272 | #dd.xz = 250.218 |
|---|
| 273 | # |
|---|
| 274 | #dd.update() |
|---|
| 275 | #dd.printmcd() |
|---|
| 276 | # |
|---|
| 277 | ###4 35,564 199,486 60,5883 14,9949 12,2972688101353 125 895,247012611329 |
|---|
| 278 | # |
|---|
| 279 | ### CASE 4 |
|---|
| 280 | #dd.lat = 35.564 |
|---|
| 281 | #dd.lon = 199.486 |
|---|
| 282 | #dd.xdate = 60.5883 |
|---|
| 283 | #dd.loct = 14.9949 |
|---|
| 284 | #dd.xz = 895.247 |
|---|
| 285 | # |
|---|
| 286 | #dd.update() |
|---|
| 287 | #dd.printmcd() |
|---|
| 288 | # |
|---|
| 289 | ###5 68,346 234,396 142,828 15,1777 13,4079899322222 83 1128,06581216829 |
|---|
| 290 | ###6 68,316 234,462 142,828 15,1819 16,1939071396022 85 582,624074786015 |
|---|
| 291 | ###7 68,302 234,144 142,828 15,1607 17,3354121022885 112 262,299040101764 |
|---|
| 292 | # |
|---|
| 293 | ### CASE 5 |
|---|
| 294 | #dd.lat = 68.32 |
|---|
| 295 | #dd.lon = 234.3 |
|---|
| 296 | #dd.xdate = 142.828 |
|---|
| 297 | #dd.loct = 15.2 |
|---|
| 298 | # |
|---|
| 299 | #dd.xz = 1128.066 |
|---|
| 300 | #dd.update() |
|---|
| 301 | #dd.printmcd() |
|---|
| 302 | # |
|---|
| 303 | #dd.xz = 582.624 |
|---|
| 304 | #dd.update() |
|---|
| 305 | #dd.printmcd() |
|---|
| 306 | # |
|---|
| 307 | #dd.xz = 262.299 |
|---|
| 308 | #dd.update() |
|---|
| 309 | #dd.printmcd() |
|---|
| 310 | # |
|---|
| 311 | # |
|---|
| 312 | ###19 27,055 129,614 13,5818 14,444 42,6488205391137 302 1395,96076100437 |
|---|
| 313 | # |
|---|
| 314 | ### CASE 19 |
|---|
| 315 | #dd.lat = 27.055 |
|---|
| 316 | #dd.lon = 129.614 |
|---|
| 317 | #dd.xdate = 13.5818 |
|---|
| 318 | #dd.loct = 14.444 |
|---|
| 319 | # |
|---|
| 320 | #dd.xz = 1395.961 |
|---|
| 321 | #dd.update() |
|---|
| 322 | #dd.printmcd() |
|---|
| 323 | # |
|---|
| 324 | #dd.lats = 20. |
|---|
| 325 | #dd.late = 30. |
|---|
| 326 | #dd.lons = 125. |
|---|
| 327 | #dd.lone = 135. |
|---|
| 328 | #dd.map2d(["wind"],incwind=True,fixedlt=True) |
|---|
| 329 | #mpl.savefig("case19.png",dpi=110,bbox_inches='tight',pad_inches=0.4) |
|---|
| 330 | # |
|---|
| 331 | # |
|---|
