| 1 | #! /usr/bin/env python |
|---|
| 2 | from ppclass import pp,ncattr |
|---|
| 3 | from ppplot import plot2d,save,writeascii |
|---|
| 4 | import numpy as np |
|---|
| 5 | from scipy.ndimage.measurements import minimum_position |
|---|
| 6 | import matplotlib.pyplot as mpl |
|---|
| 7 | import ppcompute |
|---|
| 8 | |
|---|
| 9 | ## method 3 |
|---|
| 10 | #from skimage import filter,transform,feature |
|---|
| 11 | #import matplotlib.patches as mpatches |
|---|
| 12 | ########### |
|---|
| 13 | |
|---|
| 14 | ############################################################################### |
|---|
| 15 | # FIND DUST DEVILS |
|---|
| 16 | # filefile --> file |
|---|
| 17 | # dx --> horizontal resolution |
|---|
| 18 | # dt_out --> frequency of outputs (s) |
|---|
| 19 | # lt_start --> starting local time |
|---|
| 20 | # halolim --> limit for halos |
|---|
| 21 | # method --> method used to detect vortices |
|---|
| 22 | # plotplot --> plot or not |
|---|
| 23 | # save --> save results or not |
|---|
| 24 | ############################################################################### |
|---|
| 25 | def finddd(filefile,\ |
|---|
| 26 | timelist=None,\ |
|---|
| 27 | dt_out=50.,\ |
|---|
| 28 | lt_start=8.,\ |
|---|
| 29 | halolim=None,\ |
|---|
| 30 | method=1,\ |
|---|
| 31 | plotplot=False,\ |
|---|
| 32 | save=True): |
|---|
| 33 | |
|---|
| 34 | ############################################################################### |
|---|
| 35 | ########################## FOR METHOD 1 FOR METHOD 1 ########################## |
|---|
| 36 | ## FACLIST : how many sigmas below mean we start to consider this could be a vortex |
|---|
| 37 | ## ... < 3.0 dubious low-intensity minima are caught |
|---|
| 38 | ## ... 3 ~ 3.1 is a bit too low but helps? because size of pressure drop could be an underestimate of actual dust devil |
|---|
| 39 | ## ... 3.2 ~ 3.3 is probably right, this is the one we choose for exploration [3.25] |
|---|
| 40 | ## ... NB: values 3.2 to 3.6 yields the same number of devils but measured sizes could vary a bit |
|---|
| 41 | ## ... > 3.6 makes strong minima disappear... especially >3.8 |
|---|
| 42 | ## ... EVENTUALLY 3.5-4 is better for problematic cases |
|---|
| 43 | ############################################################################### |
|---|
| 44 | faclist = [3.75] |
|---|
| 45 | |
|---|
| 46 | ## (see below) NEIGHBOR_FAC is the multiple of std used to evaluate size |
|---|
| 47 | ## --> 1: limit not discriminative enough. plus does not separate neighbouring vortices. |
|---|
| 48 | ## ... but interesting: gives an exponential law (because vortices are artificially merged?) |
|---|
| 49 | ## --> 2.7: very good for method 1. corresponds usually to ~0.3 |
|---|
| 50 | ## --> 2: so-so. do not know what to think. |
|---|
| 51 | neighbor_fac = 2.7 |
|---|
| 52 | |
|---|
| 53 | ############################################################################### |
|---|
| 54 | if save: |
|---|
| 55 | myfile1 = open(filefile+'m'+str(method)+'_'+'1.txt', 'w') |
|---|
| 56 | myfile2 = open(filefile+'m'+str(method)+'_'+'2.txt', 'w') |
|---|
| 57 | ############################################################################### |
|---|
| 58 | |
|---|
| 59 | ## get the resolution within the file |
|---|
| 60 | dx = ncattr(filefile,'DX') ; print "resolution in meters is: ",dx |
|---|
| 61 | ## if no halolim is given, guess it from resolution |
|---|
| 62 | if halolim is None: |
|---|
| 63 | extentlim = 2000. # the putative maximum extent of a vortex in m |
|---|
| 64 | halolim = extentlim / dx |
|---|
| 65 | print "maximum halo size is: ",halolim |
|---|
| 66 | |
|---|
| 67 | ## mean and std calculations |
|---|
| 68 | ## -- std is used in both methods for limits |
|---|
| 69 | ## -- mean is only used in method 1 |
|---|
| 70 | print "calculate mean and std, please wait." |
|---|
| 71 | ## -- get time series of 2D surface pressure |
|---|
| 72 | psfc = pp(file=filefile,var="PSFC").getf() |
|---|
| 73 | ## -- calculate mean and standard deviation |
|---|
| 74 | ## -- ... calculating std at all time is not right! |
|---|
| 75 | ## -- ... for mean value though, similar results with both methods |
|---|
| 76 | mean = np.mean(psfc,dtype=np.float64) |
|---|
| 77 | std = np.std(psfc,dtype=np.float64) |
|---|
| 78 | damax = np.max(psfc) |
|---|
| 79 | damin = np.min(psfc) |
|---|
| 80 | ## some information about inferred limits |
|---|
| 81 | print "**************************************************************" |
|---|
| 82 | print "MEAN",mean |
|---|
| 83 | print "STD",std |
|---|
| 84 | print "LIMIT FOR PRESSURE MINIMUM:",-np.array(faclist)*std |
|---|
| 85 | print "LIMIT FOR EVALUATING SIZE OF A GIVEN LOCAL MINIMUM",-neighbor_fac*std |
|---|
| 86 | print "**************************************************************" |
|---|
| 87 | |
|---|
| 88 | # if no timelist is given, take them all |
|---|
| 89 | if timelist is None: |
|---|
| 90 | sizet = psfc.shape[0] |
|---|
| 91 | print "treat all time values: ",sizet |
|---|
| 92 | timelist = range(0,sizet-1,1) |
|---|
| 93 | |
|---|
| 94 | ## LOOP ON TIME |
|---|
| 95 | for time in timelist: |
|---|
| 96 | |
|---|
| 97 | ## get 2D surface pressure at a given time |
|---|
| 98 | ## (this is actually so quick we don't use psfc above) |
|---|
| 99 | psfc2d = pp(file=filefile,var="PSFC",t=time).getf() |
|---|
| 100 | #ustm = pp(file=filefile,var="USTM",t=time).getf() |
|---|
| 101 | |
|---|
| 102 | ## MAIN ANALYSIS. LOOP ON FAC. OR METHOD. |
|---|
| 103 | for fac in faclist: |
|---|
| 104 | #fac = 3.75 |
|---|
| 105 | #for method in [2,1]: |
|---|
| 106 | |
|---|
| 107 | ## initialize arrays |
|---|
| 108 | tabij = [] ; tabsize = [] ; tabdrop = [] |
|---|
| 109 | tabijcenter = [] ; tabijvortex = [] ; tabijnotconv = [] ; tabdim = [] |
|---|
| 110 | |
|---|
| 111 | ################ FIND RELEVANT POINTS TO BE ANALYZED |
|---|
| 112 | ## lab is 1 for points to be treated by minimum_position routine |
|---|
| 113 | ## we set elements at 1 where pressure is under mean-fac*std |
|---|
| 114 | ## otherwise we set to 0 because this means we are close enough to mean pressure (background) |
|---|
| 115 | lab = np.zeros(psfc2d.shape) |
|---|
| 116 | if method == 1: |
|---|
| 117 | # method 1: standard deviation |
|---|
| 118 | lab[np.where(psfc2d < mean-fac*std)] = 1 |
|---|
| 119 | elif method == 2: |
|---|
| 120 | # method 2: polynomial fit |
|---|
| 121 | # ... tried smooth but too difficult, not accurate and too expensive |
|---|
| 122 | deg = 5 #plutot bien (deg 10 ~pareil) mais loupe les gros (~conv cell) |
|---|
| 123 | #deg = 2 #pas mal mais false positive (pareil 3-4 mm si un peu mieux) |
|---|
| 124 | # #(OK now with fixing the false positive bug) |
|---|
| 125 | nx = psfc2d.shape[1] ; ny = psfc2d.shape[0] |
|---|
| 126 | xxx = np.array(range(nx)) ; yyy = np.array(range(ny)) |
|---|
| 127 | anopsfc2d = psfc2d*0. ; polypsfc2d = psfc2d*0. |
|---|
| 128 | for iii in range(0,nx,1): |
|---|
| 129 | poly = np.poly1d(np.polyfit(yyy,psfc2d[iii,:],deg)) |
|---|
| 130 | polypsfc2d[iii,:] = poly(yyy) |
|---|
| 131 | for jjj in range(0,ny,1): |
|---|
| 132 | poly = np.poly1d(np.polyfit(xxx,psfc2d[:,jjj],deg)) |
|---|
| 133 | polypsfc2d[:,jjj] = 0.5*polypsfc2d[:,jjj] + 0.5*poly(xxx) |
|---|
| 134 | ## smooth a little to avoid 'crosses' (plus, this removes PBC problems) |
|---|
| 135 | polypsfc2d = ppcompute.smooth2diter(polypsfc2d,n=deg) |
|---|
| 136 | # compute anomaly and find points to be explored |
|---|
| 137 | anopsfc2d = psfc2d - polypsfc2d |
|---|
| 138 | limlim = fac*std ## same as method 1 |
|---|
| 139 | lab[np.where(anopsfc2d < -limlim)] = 1 |
|---|
| 140 | elif method == 3: |
|---|
| 141 | # method 3 : find centers of circle features using a Hough transform |
|---|
| 142 | # initialize the array containing point to be further analyzed |
|---|
| 143 | lab = np.zeros(psfc2d.shape) |
|---|
| 144 | ## prepare the field to be analyzed by the Hough transform |
|---|
| 145 | ## normalize it in an interval [-1,1] |
|---|
| 146 | field = psfc2d |
|---|
| 147 | # ... test 1. local max / min used for normalization. |
|---|
| 148 | mmax = np.max(field) ; mmin = np.min(field) ; dasigma = 2.5 |
|---|
| 149 | ## ... test 2. global max / min used for normalization. |
|---|
| 150 | #mmax = damax ; mmin = damin ; dasigma = 1.0 #1.5 trop restrictif |
|---|
| 151 | spec = 2.*((field-mmin)/(mmax-mmin) - 0.5) |
|---|
| 152 | # perform an edge detection on the field |
|---|
| 153 | # ... returns an array with True on edges and False outside |
|---|
| 154 | # http://sciunto.wordpress.com/2013/03/01/detection-de-cercles-par-une-transformation-de-hough-dans-scikit-image/ |
|---|
| 155 | edges = filter.canny(filter.sobel(spec),sigma=dasigma) |
|---|
| 156 | # initialize plot for checks |
|---|
| 157 | if plotplot: |
|---|
| 158 | fig, ax = mpl.subplots(ncols=1, nrows=1, figsize=(10,10)) |
|---|
| 159 | ax.imshow(field, cmap=mpl.cm.gray) |
|---|
| 160 | # detect circle with radius 3dx. works well. |
|---|
| 161 | # use an Hough circle transform |
|---|
| 162 | radii = np.array([3]) |
|---|
| 163 | hough_res = transform.hough_circle(edges, radii) |
|---|
| 164 | # analyze results of the Hough transform |
|---|
| 165 | nnn = 0 |
|---|
| 166 | sigselec = neighbor_fac |
|---|
| 167 | #sigselec = 3. |
|---|
| 168 | for radius, h in zip(radii, hough_res): |
|---|
| 169 | # number of circle features to keep |
|---|
| 170 | # ... quite large. but we want to be sure not to miss anything. |
|---|
| 171 | nup = 30 |
|---|
| 172 | maxima = feature.peak_local_max(h, num_peaks=nup) |
|---|
| 173 | # loop on detected circle features |
|---|
| 174 | for maximum in maxima: |
|---|
| 175 | center_x, center_y = maximum - radii.max() |
|---|
| 176 | # nup is quite high so there are false positives. |
|---|
| 177 | # ... but those are easy to detect |
|---|
| 178 | # ... if pressure drop is unclear (or inexistent) |
|---|
| 179 | # ... we do not take the point into account for further analysis |
|---|
| 180 | # ... NB: for inspection give red vs. green color to displayed circles |
|---|
| 181 | diag = field[center_x,center_y] - (mean-sigselec*std) |
|---|
| 182 | if diag < 0: |
|---|
| 183 | col = 'green' |
|---|
| 184 | nnn = nnn + 1 |
|---|
| 185 | lab[center_x,center_y] = 1 |
|---|
| 186 | else: |
|---|
| 187 | col = 'red' |
|---|
| 188 | # draw circles |
|---|
| 189 | if plotplot: |
|---|
| 190 | circ = mpatches.Circle((center_y, center_x), radius,fill=False, edgecolor=col, linewidth=2) |
|---|
| 191 | ax.add_patch(circ) |
|---|
| 192 | if plotplot: |
|---|
| 193 | mpl.title(str(nnn)+" vortices") ; mpl.show() |
|---|
| 194 | mpl.close() |
|---|
| 195 | |
|---|
| 196 | ## while there are still points to be analyzed... |
|---|
| 197 | while 1 in lab: |
|---|
| 198 | ## ... get the point with the minimum field values |
|---|
| 199 | if method == 1 or method == 3: |
|---|
| 200 | ij = minimum_position(psfc2d,labels=lab) |
|---|
| 201 | elif method == 2: |
|---|
| 202 | ij = minimum_position(anopsfc2d,labels=lab) |
|---|
| 203 | ## ... store the indexes of the point in tabij |
|---|
| 204 | tabij.append(ij) |
|---|
| 205 | ## ... remove the point from labels to be further explored by minimum_position |
|---|
| 206 | lab[ij] = 0 |
|---|
| 207 | |
|---|
| 208 | ################ GET SIZES BASED ON THOSE FOUND POINTS |
|---|
| 209 | ## reslab is the same as lab |
|---|
| 210 | ## except for scanning purpose we keep the information |
|---|
| 211 | ## about how went the detection |
|---|
| 212 | ## --> and we set a lower fac |
|---|
| 213 | ## --> above a high fac is good to catch only strong vortices |
|---|
| 214 | ## --> but here a casual fac=3 is better to get accurate sizes |
|---|
| 215 | ## --> or even lower as shown by plotting reslab |
|---|
| 216 | reslab = np.zeros(psfc2d.shape) |
|---|
| 217 | if method == 1 or method == 3: |
|---|
| 218 | reslab[np.where(psfc2d < mean-neighbor_fac*std)] = 1 |
|---|
| 219 | elif method == 2: |
|---|
| 220 | reslab[np.where(anopsfc2d < -neighbor_fac*std)] = 1 |
|---|
| 221 | |
|---|
| 222 | ## initialize halomax and while loop |
|---|
| 223 | ## HALOMAX : maximum halo defined around a minima to evaluate the size |
|---|
| 224 | ## ... halomax must be large enough to encompass a vortex |
|---|
| 225 | ## ... but not too large otherwise neighboring vortex are caught |
|---|
| 226 | halomax = 3 ; notconv = 9999 ; yorgl = 9999 |
|---|
| 227 | |
|---|
| 228 | ## WHILE LOOP on HALOMAX exploration |
|---|
| 229 | while ( notconv > 0 and halomax < halolim and yorgl != 0 ): |
|---|
| 230 | # now browse through all points caught in previous loop |
|---|
| 231 | for ij in tabij: |
|---|
| 232 | ## ... OK. take each indexes found before with minimum_position |
|---|
| 233 | i,j = ij[0],ij[1] |
|---|
| 234 | ## EITHER |
|---|
| 235 | ## ... if reslab is already 0, we do not have to do anything |
|---|
| 236 | ## ... because this means point is already part of detected vortex |
|---|
| 237 | ## OR |
|---|
| 238 | ## ... if the ij couple is already in a previously detected vortex |
|---|
| 239 | ## ... we don't actually need to consider it again |
|---|
| 240 | ## ... this is necessary otherwise (sometimes a lot) of false positives |
|---|
| 241 | if reslab[i,j] <= 0 or ij in tabijvortex: |
|---|
| 242 | pass |
|---|
| 243 | else: |
|---|
| 244 | ## ... then define a growing halo around this point |
|---|
| 245 | ## ... we make the halo grow until convergence of |
|---|
| 246 | ## ... the number of under-limit points (i.e. with lab=1) within halo |
|---|
| 247 | ## ... which means the whole vortex is encompassed |
|---|
| 248 | ## ... we start with a halo of 1 point around minimum point |
|---|
| 249 | ## ... and we end when we reached halomax which likely means dubious case |
|---|
| 250 | halo = 1 ; nmesh = -9999. ; prevmesh = 9999. ; notconverged = False |
|---|
| 251 | while nmesh != prevmesh and halo <= halomax: |
|---|
| 252 | ## ... save the number of vortex points calculated at previous iteration |
|---|
| 253 | prevmesh = nmesh |
|---|
| 254 | ## ... define a halo around the minimum point |
|---|
| 255 | minx,maxx,miny,maxy = i-halo,i+halo+1,j-halo,j+halo+1 |
|---|
| 256 | ## ... treat the boundary case (TBD: periodic boundary conditions) |
|---|
| 257 | if minx < 0: minx = 0 |
|---|
| 258 | if miny < 0: miny = 0 |
|---|
| 259 | if maxx > psfc2d.shape[0]: maxx = psfc2d.shape[0] |
|---|
| 260 | if maxy > psfc2d.shape[1]: maxy = psfc2d.shape[1] |
|---|
| 261 | ## ... define the patch, made of the halo of points |
|---|
| 262 | patch = reslab[minx:maxx,miny:maxy] |
|---|
| 263 | ## ... count how many 1 are inside this patch |
|---|
| 264 | ## ... these are points for which value is >0 |
|---|
| 265 | ## ... because close to mean is 0 and already caught is <0 |
|---|
| 266 | ## ... and not converged are >1 so still to be scanned |
|---|
| 267 | nmesh = len(np.where(patch >= 1)[0]) |
|---|
| 268 | ## ... in case widening halo adds 1-2 points only |
|---|
| 269 | ## ... we consider vortex is encompassed |
|---|
| 270 | ## ... this helps to separate neighboring vortices |
|---|
| 271 | ## ... and only yields marginal inaccuracies |
|---|
| 272 | if halo > 3 and abs(nmesh-prevmesh)<=2: prevmesh = nmesh |
|---|
| 273 | ## ... if with halo=1 we caught a 1-point vortex, end here because spurious |
|---|
| 274 | if halo == 1 and nmesh == 1: prevmesh = nmesh |
|---|
| 275 | ## ... with the last halo=halomax test, we can detect not-converged cases |
|---|
| 276 | if halo == halomax and nmesh != prevmesh: notconverged = True |
|---|
| 277 | ## ... increment halo size before the next iteration |
|---|
| 278 | halo = halo + 1 |
|---|
| 279 | ## ... now we got the grid points encompassed by the vortex in nmesh |
|---|
| 280 | ## ... then to continue scanning we store results in reslab |
|---|
| 281 | if notconverged: |
|---|
| 282 | ## multiply reslab by 2. |
|---|
| 283 | ## --> if it is >1, point could be part of another vortex found later |
|---|
| 284 | ## --> if it is <0, point is already within a caught vortex, so reslab should remain <0 |
|---|
| 285 | ## --> if it is =0, this should remain 0 because close to average |
|---|
| 286 | reslab[i,j] = reslab[i,j]*2 |
|---|
| 287 | tabijnotconv.append(ij) |
|---|
| 288 | else: |
|---|
| 289 | ## OK. this is most likely an actual vortex. we get the drop. |
|---|
| 290 | ## we multiply by mesh area, then square to get approx. size of vortex |
|---|
| 291 | if method == 1 or method ==3: |
|---|
| 292 | drop = -psfc2d[i,j]+mean |
|---|
| 293 | else: |
|---|
| 294 | drop = -anopsfc2d[i,j] |
|---|
| 295 | size = int(np.sqrt(nmesh*dx*dx)) |
|---|
| 296 | ## if vortex is too small (i.e. too close to mesh grid resolution) we don't store information |
|---|
| 297 | ## ... we just remove the patch from labels to be further explored |
|---|
| 298 | facdx = 2. |
|---|
| 299 | if size < facdx*dx: |
|---|
| 300 | reslab[minx:maxx,miny:maxy] = 0 |
|---|
| 301 | ## otherwise it is a VORTEX! we store info in arrays |
|---|
| 302 | else: |
|---|
| 303 | ## we put the vortex points into tabijvortex so that |
|---|
| 304 | ## a (i,j) couple into a vortex is not considered in further explorations |
|---|
| 305 | ## -- also we evaluate the x-size and y-size of the vortex (max distance to center) |
|---|
| 306 | ## -- this can be useful to detect not-so-round fake vortices (convective gusts?) |
|---|
| 307 | maxw = -9999 ; maxh = -9999 ; maxu = -9999 |
|---|
| 308 | for iii in range(minx,maxx): |
|---|
| 309 | for jjj in range(miny,maxy): |
|---|
| 310 | if reslab[iii,jjj] != 0: |
|---|
| 311 | tabijvortex.append((iii,jjj)) |
|---|
| 312 | width = np.abs(iii-ij[0]) |
|---|
| 313 | if width > maxw: maxw = width |
|---|
| 314 | height = np.abs(jjj-ij[1]) |
|---|
| 315 | if height > maxh: maxh = height |
|---|
| 316 | #ustar = ustm[iii,jjj] |
|---|
| 317 | #if ustar > maxu: maxu = ustar |
|---|
| 318 | #print size,drop,ustar |
|---|
| 319 | ## store info in dedicated arrays |
|---|
| 320 | tabdim.append((maxw*dx,maxh*dx)) |
|---|
| 321 | tabsize.append(size) |
|---|
| 322 | tabdrop.append(drop) |
|---|
| 323 | tabijcenter.append(ij) |
|---|
| 324 | #print "... VORTEX!!!! size %.0f drop %.1f coord %.0f %.0f" % (size,drop,i,j) |
|---|
| 325 | ## we remove the patch from labels to be further explored |
|---|
| 326 | ## we multiply reslab by -1 to plot detection maps |
|---|
| 327 | reslab[minx:maxx,miny:maxy] = patch*-1 |
|---|
| 328 | |
|---|
| 329 | ## count how many points are not converged and left to be analyzed |
|---|
| 330 | notconv = len(np.where(reslab > 1)[0]) |
|---|
| 331 | yorgl = len(np.where(reslab == 1)[0]) |
|---|
| 332 | |
|---|
| 333 | ## increment halomax |
|---|
| 334 | ## to speed-up the increment is slightly increasing with considered halomax |
|---|
| 335 | halomax = halomax + halomax / 2 |
|---|
| 336 | |
|---|
| 337 | ## just for simpler plots. |
|---|
| 338 | reslab[reslab > 2] = 4 |
|---|
| 339 | reslab[reslab < -2] = -4 |
|---|
| 340 | |
|---|
| 341 | ## give some info to the user |
|---|
| 342 | if len(tabsize) > 0: |
|---|
| 343 | nvortex = len(tabsize) |
|---|
| 344 | maxsize = np.max(tabsize) |
|---|
| 345 | maxdrop = np.max(tabdrop) |
|---|
| 346 | else: |
|---|
| 347 | nvortex = 0 |
|---|
| 348 | maxsize = 0 |
|---|
| 349 | maxdrop = 0. |
|---|
| 350 | notconv = len(np.where(reslab > 1)[0]) |
|---|
| 351 | print "t=%3.0f / n=%2.0f / s_max=%4.0f / d_max=%4.1f / halo_out=%3.0f / notconvp=%3.1f" \ |
|---|
| 352 | % (time,nvortex,maxsize,maxdrop,halomax,100.*notconv/float(reslab.size)) |
|---|
| 353 | |
|---|
| 354 | ## save results in a text file |
|---|
| 355 | if save: |
|---|
| 356 | # convert t in local time |
|---|
| 357 | ttt = lt_start + time*dt_out/3700. |
|---|
| 358 | # write files |
|---|
| 359 | myfile2.write( "%5.2f ; %5.0f ; %5.0f ; %7.2f\n" % (ttt,nvortex,maxsize,maxdrop) ) |
|---|
| 360 | for iii in range(len(tabsize)): |
|---|
| 361 | myfile1.write( "%5.2f ; %5.0f ; %7.2f ; %5.0f ; %5.0f\n" \ |
|---|
| 362 | % (ttt,tabsize[iii],tabdrop[iii],tabdim[iii][0],tabdim[iii][1]) ) |
|---|
| 363 | |
|---|
| 364 | #### PLOT PLOT PLOT PLOT |
|---|
| 365 | if (nvortex>0 and plotplot) or (nvortex>0 and maxsize > 800.): |
|---|
| 366 | mpl.figure(figsize=(12,8)) |
|---|
| 367 | myplot = plot2d() |
|---|
| 368 | myplot.x = np.array(range(psfc2d.shape[1]))*dx/1000. |
|---|
| 369 | myplot.y = np.array(range(psfc2d.shape[0]))*dx/1000. |
|---|
| 370 | myplot.title = str(nvortex)+" vortices found (indicated diameter / pressure drop)" |
|---|
| 371 | myplot.xlabel = "x distance (km)" |
|---|
| 372 | myplot.ylabel = "y distance (km)" |
|---|
| 373 | if method > 0: |
|---|
| 374 | #if method == 1: |
|---|
| 375 | #myplot.field = ustm |
|---|
| 376 | myplot.f = psfc2d |
|---|
| 377 | #myplot.field = polypsfc2d |
|---|
| 378 | myplot.vmin = -2.*std |
|---|
| 379 | myplot.vmax = +2.*std |
|---|
| 380 | myplot.vmin = mean - 6.*std |
|---|
| 381 | myplot.vmax = mean + 6.*std |
|---|
| 382 | else: |
|---|
| 383 | myplot.field = anopsfc2d |
|---|
| 384 | myplot.vmin = -1.5 |
|---|
| 385 | myplot.vmax = 0.5 |
|---|
| 386 | myplot.fmt = "%.1f" |
|---|
| 387 | myplot.div = 20 |
|---|
| 388 | myplot.colorbar = "spectral" |
|---|
| 389 | myplot.make() |
|---|
| 390 | |
|---|
| 391 | ### ANNOTATIONS |
|---|
| 392 | for iii in range(len(tabsize)): |
|---|
| 393 | ij = tabijcenter[iii] |
|---|
| 394 | coord1 = ij[1]*dx/1000. |
|---|
| 395 | coord2 = ij[0]*dx/1000. |
|---|
| 396 | txt = "%.0f/%.1f" % (tabsize[iii],tabdrop[iii]) |
|---|
| 397 | mpl.annotate(txt,xy=(coord1,coord2), |
|---|
| 398 | xytext=(-10,-30),textcoords='offset points',ha='center',va='bottom',\ |
|---|
| 399 | bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),\ |
|---|
| 400 | arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3',color='red'),\ |
|---|
| 401 | size='small') |
|---|
| 402 | |
|---|
| 403 | ###show detection |
|---|
| 404 | #lev = [-4,-2,-1,0,1,2,4] |
|---|
| 405 | #lev = [-1,0] |
|---|
| 406 | #mpl.contourf(myplot.absc,myplot.ordi,reslab,alpha=0.9,cmap=mpl.cm.get_cmap("binary_r"),levels=lev) |
|---|
| 407 | |
|---|
| 408 | ### SHOW OR SAVE IN FILE |
|---|
| 409 | mpl.show() |
|---|
| 410 | #save(mode="png",filename="detectm"+"_"+str(time)+"_"+str(method),folder="detect/",includedate=False) |
|---|
| 411 | |
|---|
| 412 | ## close data files |
|---|
| 413 | if save: |
|---|
| 414 | myfile1.close() |
|---|
| 415 | myfile2.close() |
|---|
