[2] | 1 | ;------------------------------------------------------------ |
---|
| 2 | ;------------------------------------------------------------ |
---|
| 3 | ;------------------------------------------------------------ |
---|
| 4 | ;+ |
---|
| 5 | ; NAME: |
---|
| 6 | ; |
---|
| 7 | ; PURPOSE: |
---|
| 8 | ; |
---|
| 9 | ; CATEGORY: |
---|
| 10 | ; |
---|
| 11 | ; CALLING SEQUENCE: |
---|
| 12 | ; |
---|
| 13 | ; INPUTS: |
---|
| 14 | ; |
---|
| 15 | ; KEYWORD PARAMETERS: |
---|
| 16 | ; |
---|
| 17 | ; OUTPUTS: |
---|
| 18 | ; |
---|
| 19 | ; COMMON BLOCKS:common.pro |
---|
| 20 | ; |
---|
| 21 | ; SIDE EFFECTS: |
---|
| 22 | ; |
---|
| 23 | ; RESTRICTIONS: |
---|
| 24 | ; |
---|
| 25 | ; EXAMPLE: |
---|
| 26 | ; |
---|
| 27 | ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) |
---|
| 28 | ; |
---|
| 29 | ;- |
---|
| 30 | ;------------------------------------------------------------ |
---|
| 31 | ;------------------------------------------------------------ |
---|
| 32 | ;------------------------------------------------------------ |
---|
| 33 | |
---|
| 34 | PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints, BOITE = boite, TYPE = type, DIREC = direc |
---|
| 35 | @common |
---|
| 36 | ; |
---|
| 37 | ;--------------------- |
---|
| 38 | array = litchamp(field) |
---|
| 39 | ;--------------------- |
---|
| 40 | ; definition de boite en fonction de endpoints |
---|
| 41 | ; puis redefinition du domaine |
---|
| 42 | boite2d = [min([endpoints[0], endpoints[2]]), max([endpoints[0], endpoints[2]]), min([endpoints[1], endpoints[3]]), max([endpoints[1], endpoints[3]])] |
---|
| 43 | ; |
---|
| 44 | minprof = 0 |
---|
| 45 | profdefault = 200 |
---|
| 46 | if n_elements(type) EQ 0 then type = 'nothing' |
---|
| 47 | Case N_Elements(Boite) OF |
---|
| 48 | 1:boite=[boite2d, minprof, boite[0]] |
---|
| 49 | 2:boite=[boite2d, boite[0],boite[1]] |
---|
| 50 | Else:$ |
---|
| 51 | if strpos(type, 'z') NE -1 THEN boite=[boite2d, minprof, profdefault] $ |
---|
| 52 | ELSE boite=[boite2d, prof1, prof2] |
---|
| 53 | ENDCASE |
---|
| 54 | if strpos(type, 'z') NE -1 THEN BEGIN |
---|
| 55 | !y.range = [boite[5], boite[4]] ;on garde les yranges (axe z) avant de changer la boite. |
---|
| 56 | profmax = boite[5] |
---|
| 57 | vraiboite = boite |
---|
| 58 | if vargrid EQ 'W' then gdep = gdepw ELSE gdep = gdept |
---|
| 59 | ; check with vertical grid limits (nearest level) |
---|
| 60 | gwork = gdep |
---|
| 61 | ; check the increse or decrese of the z axis |
---|
| 62 | IF gwork[1] LE gwork[0] THEN gwork = reverse(gdep, 1) |
---|
| 63 | niveauprof = where(gwork ge boite[5]) & niveauprof = niveauprof[0] |
---|
| 64 | if niveauprof NE -1 then boite[5] = gwork[niveauprof]+1 |
---|
| 65 | ENDIF |
---|
| 66 | domdef, boite, GRILLE=[vargrid], /findalways, _extra = ex |
---|
| 67 | ; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef |
---|
| 68 | lon1 = min([endpoints[0], endpoints[2]]) |
---|
| 69 | lon2 = max([endpoints[0], endpoints[2]]) |
---|
| 70 | lat1 = min([endpoints[1], endpoints[3]]) |
---|
| 71 | lat2 = max([endpoints[1], endpoints[3]]) |
---|
| 72 | ; on recupere la grille sur la boite |
---|
| 73 | grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz |
---|
| 74 | ;--------------------- |
---|
| 75 | ; on definit la triangulation qui va nous permetre de determiner la |
---|
| 76 | ; section. on la recalcule car elle doit etre definie sur la terre |
---|
| 77 | ; aussi bien que sur la mer. suivant le sens de la section -plutot |
---|
| 78 | ; longitude ou plutot latitude- on definit la facon de trianguler |
---|
| 79 | if strpos(type, 'x') NE -1 then BEGIN |
---|
| 80 | downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] |
---|
| 81 | tri = definetri(nx, ny, (downward)[*]) |
---|
| 82 | ENDIF ELSE tri = definetri(nx, ny) |
---|
| 83 | ;--------------------- |
---|
| 84 | ; equation de la droite suivant laquelle on fait la section |
---|
| 85 | ;--------------------- |
---|
| 86 | abc = linearequation(endpoints[0:1], endpoints[2:3]) |
---|
| 87 | glamtri = glam[tri] |
---|
| 88 | gphitri = gphi[tri] |
---|
| 89 | ; quels sont les points de la triangulation qui sont au dessus et au |
---|
| 90 | ; dessous de la droite ? |
---|
| 91 | if abc[1] NE 0 THEN $ |
---|
| 92 | test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $ |
---|
| 93 | ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) |
---|
| 94 | |
---|
| 95 | zero123 = total(test, 1) |
---|
| 96 | ; to keep: triangles de la triangulation qui sont a cheval sur la droite. |
---|
| 97 | tokeep1 = where(zero123 EQ 1) |
---|
| 98 | tokeep2 = where(temporary(zero123) EQ 2) |
---|
| 99 | tokeep = [tokeep1, tokeep2] |
---|
| 100 | |
---|
| 101 | test = test[*, tokeep] |
---|
| 102 | tri = tri[*, tokeep] |
---|
| 103 | ; quel est le sommet du triangle qui est seul d''un cote de la droite? |
---|
| 104 | single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) |
---|
| 105 | single1 = single1-(single1/3)*3 |
---|
| 106 | single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0) |
---|
| 107 | single2 = single2-(single2/3)*3 |
---|
| 108 | |
---|
| 109 | undefine, tokeep |
---|
| 110 | undefine, tokeep1 |
---|
| 111 | undefine, tokeep2 |
---|
| 112 | undefine, test |
---|
| 113 | |
---|
| 114 | single = [temporary(single1), temporary(single2)] |
---|
| 115 | ; points1 le point du triangles qui est seul d''un cote de la droite. |
---|
| 116 | ; point2 l''autre point du triangle de l''autre cote de la droite |
---|
| 117 | point1 = [single, single] |
---|
| 118 | point2 = [single EQ 0, 1 + (single LE 1)] |
---|
| 119 | |
---|
| 120 | undefine, single |
---|
| 121 | |
---|
| 122 | ntri = (size(tri))[2] |
---|
| 123 | index = [lindgen(ntri), lindgen(ntri)] |
---|
| 124 | |
---|
| 125 | points1 = tri[point1, index] |
---|
| 126 | points2 = tri[point2, temporary(index)] |
---|
| 127 | ; points : complexe contenant les couples de points de part et |
---|
| 128 | ; d''autre de la droite. Ils faut supprimer les doublons. |
---|
| 129 | points = dcomplex(points1, points2) |
---|
| 130 | points = points[uniq(points, sort(points))] |
---|
| 131 | symetrique = dcomplex(imaginary(points), double(points)) |
---|
| 132 | points = points[where(points-shift(temporary(symetrique), 1) NE 0)] |
---|
| 133 | ; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite. |
---|
| 134 | ; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite |
---|
| 135 | points1 = complex(glam[ double(points)], gphi[ double(points)]) |
---|
| 136 | points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) |
---|
| 137 | ; droites les equations des droites dont on cherche l''intersection |
---|
| 138 | ; avec la section. |
---|
| 139 | droites = linearequation(points1, points2) |
---|
| 140 | inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) |
---|
| 141 | |
---|
| 142 | ; les ccordonnes geographiques des points que l''on cherche sur la section. |
---|
| 143 | glamaxe = float(inter) |
---|
| 144 | gphiaxe = imaginary(inter) |
---|
| 145 | ; on les range ds l''ordre croissant entre les bornes de la section |
---|
| 146 | if strpos(type, 'x') NE -1 then BEGIN |
---|
| 147 | sort = sort(glamaxe) |
---|
| 148 | glamaxe = glamaxe[sort] |
---|
| 149 | inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) |
---|
| 150 | glamaxe = glamaxe[inbox] |
---|
| 151 | sort = sort[inbox] |
---|
| 152 | gphiaxe = gphiaxe[sort] |
---|
| 153 | ENDIF ELSE BEGIN |
---|
| 154 | sort = sort(gphiaxe) |
---|
| 155 | gphiaxe = gphiaxe[sort] |
---|
| 156 | inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) |
---|
| 157 | gphiaxe = gphiaxe[inbox] |
---|
| 158 | sort = sort[inbox] |
---|
| 159 | glamaxe = glamaxe[sort] |
---|
| 160 | ENDELSE |
---|
| 161 | points = points[sort] |
---|
| 162 | points1 = points1[sort] |
---|
| 163 | points2 = points2[sort] |
---|
| 164 | inter = inter[sort] |
---|
| 165 | poids = abs(points2-inter)/abs(points2-points1) |
---|
| 166 | |
---|
| 167 | ;--------------------- |
---|
| 168 | array = fitintobox(array) |
---|
| 169 | if array[0] EQ -1 THEN BEGIN |
---|
| 170 | res = -1 |
---|
| 171 | return |
---|
| 172 | ENDIF |
---|
| 173 | if n_elements(valmask) EQ 0 THEN valmask = 1e20 |
---|
| 174 | taille=size(array) |
---|
| 175 | if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN |
---|
| 176 | direc = 't' |
---|
| 177 | array = grossemoyenne(array, 't') |
---|
| 178 | taille=size(array) |
---|
| 179 | jpt = 1 |
---|
| 180 | ENDIF |
---|
| 181 | case 1 of |
---|
| 182 | ;---------------------------------------------------------------------------- |
---|
| 183 | ;xy |
---|
| 184 | ;---------------------------------------------------------------------------- |
---|
| 185 | taille[0] EQ 2:BEGIN |
---|
| 186 | value1 = array[double(points)] |
---|
| 187 | terre = where(value1 GT valmask/10) |
---|
| 188 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
| 189 | value2 = array[imaginary(points)] |
---|
| 190 | terre = where(value2 GT valmask/10) |
---|
| 191 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
| 192 | res = poids*value1+(1-poids)*value2 |
---|
| 193 | END |
---|
| 194 | ;---------------------------------------------------------------------------- |
---|
| 195 | ;xyz |
---|
| 196 | ;---------------------------------------------------------------------------- |
---|
| 197 | taille[0] EQ 3 AND jpt EQ 1:BEGIN |
---|
| 198 | npoints = n_elements(points) |
---|
| 199 | index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) |
---|
| 200 | value1 = array[index] |
---|
| 201 | terre = where(value1 GT valmask/10) |
---|
| 202 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
| 203 | index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz) |
---|
| 204 | value2 = array[index] |
---|
| 205 | terre = where(value2 GT valmask/10) |
---|
| 206 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
| 207 | poids = poids#replicate(1, nz) |
---|
| 208 | res = poids*value1+(1-poids)*value2 |
---|
| 209 | ; moyenne suivant z ? |
---|
| 210 | if strpos(type, 'z') EQ -1 then begin |
---|
| 211 | nan = where(finite(res) EQ 0) |
---|
| 212 | if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] |
---|
| 213 | weight = replicate(1, npoints)#e3 |
---|
| 214 | if nan[0] NE -1 then weight[nan] = !values.f_nan |
---|
| 215 | totalweight = total(weight, 2, /nan) |
---|
| 216 | zero = where(totalweight EQ 0) |
---|
| 217 | if zero[0] NE -1 then totalweight[zero] = !values.f_nan |
---|
| 218 | res = total(res*weight, 2, /nan)/totalweight |
---|
| 219 | direc = 'z'+string(byte(testvar(var=toto))) |
---|
| 220 | endif |
---|
| 221 | END |
---|
| 222 | ;---------------------------------------------------------------------------- |
---|
| 223 | ;xyt |
---|
| 224 | ;---------------------------------------------------------------------------- |
---|
| 225 | taille[0] EQ 3 AND jpt NE 1:BEGIN |
---|
| 226 | npoints = n_elements(points) |
---|
| 227 | index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) |
---|
| 228 | value1 = array[index] |
---|
| 229 | terre = where(value1 GT valmask/10) |
---|
| 230 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
| 231 | index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt) |
---|
| 232 | value2 = array[index] |
---|
| 233 | terre = where(value2 GT valmask/10) |
---|
| 234 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
| 235 | poids = poids#replicate(1, jpt) |
---|
| 236 | res = poids*value1+(1-poids)*value2 |
---|
| 237 | END |
---|
| 238 | ;---------------------------------------------------------------------------- |
---|
| 239 | ;xyzt |
---|
| 240 | ;---------------------------------------------------------------------------- |
---|
| 241 | taille[0] EQ 4:BEGIN |
---|
| 242 | npoints = n_elements(points) |
---|
| 243 | index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) |
---|
| 244 | index = reform(index, npoints, nz, jpt, /over) |
---|
| 245 | value1 = array[index] |
---|
| 246 | terre = where(value1 GT valmask/10) |
---|
| 247 | if terre[0] NE -1 then value1[terre] = !values.f_nan |
---|
| 248 | index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt) |
---|
| 249 | index = reform(index, npoints, nz, jpt, /over) |
---|
| 250 | value2 = array[index] |
---|
| 251 | terre = where(value2 GT valmask/10) |
---|
| 252 | if terre[0] NE -1 then value2[terre] = !values.f_nan |
---|
| 253 | poids = poids#replicate(1, nz*jpt) |
---|
| 254 | poids = reform(poids, npoints, nz, jpt, /over) |
---|
| 255 | res = poids*value1+(1-poids)*value2 |
---|
| 256 | ; moyenne suivant z ? |
---|
| 257 | if strpos(type, 'z') EQ -1 then begin |
---|
| 258 | nan = where(finite(res) EQ 0) |
---|
| 259 | if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] |
---|
| 260 | weight = replicate(1, npoints)#e3 |
---|
| 261 | weight = weight[*]#replicate(1, jpt) |
---|
| 262 | weight = reform(weight, npoints, nz, jpt, /over) |
---|
| 263 | if nan[0] NE -1 then weight[nan] = !values.f_nan |
---|
| 264 | totalweight = total(weight, 2, /nan) |
---|
| 265 | zero = where(totalweight EQ 0) |
---|
| 266 | if zero[0] NE -1 then totalweight[zero] = !values.f_nan |
---|
| 267 | res = total(res*weight, 2, /nan)/totalweight |
---|
| 268 | direc = 'z'+string(byte(testvar(var=toto))) |
---|
| 269 | endif |
---|
| 270 | END |
---|
| 271 | endcase |
---|
| 272 | ;--------------------- |
---|
| 273 | |
---|
| 274 | terre = where(finite(res) EQ 0) |
---|
| 275 | if terre[0] NE -1 then res[terre] = valmask |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | ; plt, tmask[*,*,0],/nocouleur, /rempli, title = '', subtitle = '' |
---|
| 279 | ; !p.title='' |
---|
| 280 | ; !p.subtitle='' |
---|
| 281 | ; dessinetri |
---|
| 282 | ; plot,[endpoints[0],endpoints[2]],[endpoints[1],endpoints[3]],/noerase,color=50 |
---|
| 283 | |
---|
| 284 | ; plot,float(points1),imaginary(points1),/noerase,color=150,psym=1 |
---|
| 285 | ; plot,float(points2),imaginary(points2),/noerase,color=150,psym=1 |
---|
| 286 | ; plot,float(inter),imaginary(inter),/noerase,color=250,psym=1 |
---|
| 287 | |
---|
| 288 | ; plot,[float(points1[0])],[imaginary(points1[0])],/noerase,color=150,psym=1 |
---|
| 289 | ; plot,[float(points2[0])],[imaginary(points2[0])],/noerase,color=150,psym=1 |
---|
| 290 | ; plot,[float(inter[0])],[imaginary(inter[0])],/noerase,color=250,psym=1 |
---|
| 291 | |
---|
| 292 | |
---|
| 293 | ;--------------------- |
---|
| 294 | return |
---|
| 295 | end |
---|