Changeset 150 for trunk/SRC/ToBeReviewed/TRIANGULATION
- Timestamp:
- 08/09/06 12:12:54 (18 years ago)
- Location:
- trunk/SRC/ToBeReviewed/TRIANGULATION
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SRC/ToBeReviewed/TRIANGULATION/ciseauxtri.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME: 6 ; 7 ; PURPOSE:vire les tableaux qui ne doivent pas etre dessines grace a 2 8 ; tests: 9 ; 1) les coins du tableau doivent etre ds la fenetre 10 ; 2) les clongeurs des cotes des triangfles exprimes en 11 ; coordonnees normalisesne doivent pas depasser une certaine 12 ; longueur seuil. 13 ; 14 ; CATEGORY: 15 ; 16 ; CALLING SEQUENCE: 17 ; 18 ; INPUTS: 19 ; 20 ; KEYWORD PARAMETERS: 21 ; 22 ; OUTPUTS: 23 ; 24 ; COMMON BLOCKS: 25 ; common.pro 26 ; 27 ; SIDE EFFECTS: 28 ; 29 ; RESTRICTIONS: 30 ; 31 ; EXAMPLE: 32 ; 33 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 5 ; @file_comments 6 ; Delete arrays which do not have to be drawn thanks to 2 tests: 7 ; 1) Corners of the array must be in the window 8 ; 2) Lenghtes of side of triangles expressed in normalized 9 ; coordinates must not surpass a sill lenght. 10 ; 11 ; @categories 12 ; 13 ; @param TRIANG 14 ; 15 ; @param GLAM 16 ; 17 ; @param GPHI 18 ; 19 ; @keyword ALL 20 ; 21 ; @keyword _EXTRA 22 ; Used to pass your keywords. 23 ; 24 ; @uses 25 ; common.pro 26 ; 27 ; @history 28 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 34 29 ; 20/2/99 30 ; 31 ; @version 32 ; $Id$ 33 ; 35 34 ;- 36 35 ;------------------------------------------------------------ 37 36 ;------------------------------------------------------------ 38 37 ;------------------------------------------------------------ 39 FUNCTION ciseauxtri, triang, glam, gphi, TOUT = tout, _EXTRA = ex38 FUNCTION ciseauxtri, triang, glam, gphi, ALL = all, _EXTRA = ex 40 39 ;--------------------------------------------------------- 41 40 ; … … 52 51 OR !map.projection EQ 18) THEN return, triang 53 52 ; 54 tempsun = systime(1) ; pour key_performance53 tempsun = systime(1) ; For key_performance 55 54 ; 56 55 taille = size(glam) … … 58 57 ny = taille[2] 59 58 60 tempdeux = systime(1) ; pour key_performance =259 tempdeux = systime(1) ; For key_performance =2 61 60 z = convert_coord(glam[*],gphi[*],/data,/to_normal) 62 61 x = z[0, *] … … 66 65 print, 'temps ciseauxtri: convert_coord data to_normal', systime(1)-tempdeux 67 66 ; 68 ; attention, suivant la projection certains points x ou y peuvent 69 ; devenir NaN (cf. points deriere la terre ds une projection 70 ; orthographique) il faut dans ce cas enlever tous les triangles qui 71 ; contiennent un de ces points 67 ; Beware, following the projection, some points x or y can become NaN 68 ; (see points behind the Earth in an orthographic projection). 69 ; In this case, we have to remove all triangle which contain one of these points. 72 70 ; 73 71 if (!map.projection LE 7 AND !map.projection NE 0) $ 74 72 OR !map.projection EQ 14 OR !map.projection EQ 15 OR !map.projection EQ 18 then begin 75 tempdeux = systime(1) ; pour key_performance =273 tempdeux = systime(1) ; For key_performance =2 76 74 ; 77 75 test = (x*y)[triang] … … 88 86 seuil = 5 < (min([nx, ny])-2) 89 87 ; 90 ; maintenant on vire les triangles dont un des cotes a une taille 91 ; superieure a 1/seuil du domaine reserve au dessin. 88 ; Now we delete tringles whose one side has a size superior to 1/sill from the domain reserved for the drawing. 92 89 ; 93 90 94 91 if keyword_set(key_periodic) then begin 95 tempdeux = systime(1) ; pour key_performance =292 tempdeux = systime(1) ; For key_performance =2 96 93 ; 97 94 xtri = x[triang] … … 111 108 print, 'temps ciseauxtri: trouver les triangles trop grands', systime(1)-tempdeux 112 109 ; 113 tempdeux = systime(1) ; pour key_performance =2110 tempdeux = systime(1) ; For key_performance =2 114 111 if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return, -1 115 112 trichanged = 1b … … 118 115 endif 119 116 ; 120 ; on supprime les triangles qui sont un peu trop hors du cadre121 ; valable qd /TOUT est active117 ; We delete all triangle which are out of the valid frame when on supprime les triangles qui sont un peu trop hors du cadre 118 ; valable qd ALL is activated 122 119 ; 123 120 if keyword_set(key_irregular) then begin 124 121 125 tempdeux = systime(1) ; pour key_performance =2122 tempdeux = systime(1) ; For key_performance =2 126 123 xtri = x[triang] 127 124 test1 = xtri GE !p.position[0] … … 148 145 ENDIF 149 146 ; 150 ; dernier tri147 ; Last sort. 151 148 ; 152 149 if keyword_set(trichanged) then BEGIN 153 150 ; 154 ; il faut verifier que les triangles que l''on garde ne 155 ; forment pas une triangulation dans laquelle 2 triangles ont un 156 ; sommet en commun sans avoir de cote de commun. 157 ; 158 ; on constitue la liste des rectangles que l''on souhaite garder (on 159 ; garde un rectangle des qu''il y a un triangle dedans) 160 ; 161 ; dans definetri, on a construit les triangles tels que le premier et 162 ; le dernier sommets soient ceux de la diagonale du rectangle definit 163 ; par le maillage. 164 ; pour retouver de quel rectangle provient un triangle, on cherche le 165 ; min de l''indice suivant x et suivant y de chaque triangle. Apres on 166 ; repasse cette indissage suivant x et y en indicage suivant nx*ny 167 ; 168 tempdeux = systime(1) ; pour key_performance =2 169 ; indices en y des rectangles 151 ; We have to check that triangles we keep do not form a triangulation 152 ; in which 2 triangles have a common summit without have a common side. 153 ; 154 ; We constitute the list of rectangles we want to keep (we keep every 155 ; rectangle containing a triangle) 156 ; 157 ; In definetri, we have construct triangles just so the first and the 158 ; last summit are those of the diagonale of the rectangle defined by 159 ; the mesh size. To find from which rectangle a triangle comes, we look 160 ; for the min of the index folowing x and following y of each triangle. 161 ; Then we go by again this indexion following x and y in an indexion 162 ; following nx*ny/ 163 ; 164 tempdeux = systime(1) ; For key_performance =2 165 ; y indexes of rectangles 170 166 indytriang = (triang[0, *]/nx) < (triang[2, *]/nx) 171 ; indices en x desrectangles167 ; x indexes of rectangles 172 168 indxtriang0 = triang[0, *]-nx*(triang[0, *]/nx) 173 169 indxtriang2 = triang[2, *]-nx*(triang[2, *]/nx) 174 170 indxmin = indxtriang0 < indxtriang2 175 ; attention dans le cas ou la grille est periodique et ou on a tous176 ; les points suivant x, les triangles qui assurent la periodicite en x177 ; ont des indices suivant x qui sont 0 et nx-1. Ils appatienent au178 ; rectangles de la colonne nx-1 et non de la colonne 0171 ; Beware in the case where the grid is periodic and where we have all points 172 ; following x, triangles which assure the periodicity in x have indexes 173 ; following x which are 0 and nx-1. They belong to rectangles of the column 174 ; nx-1 and not to column 0. 179 175 if keyword_set(key_periodic) AND nx EQ jpi then BEGIN 180 176 indxmax = indxtriang0 > indxtriang2 … … 185 181 print, 'temps ciseauxtri: liste des rectangles', systime(1)-tempdeux 186 182 ; 187 ; maintenant qu''on a cette liste on va s''assuter que l''on a pas de188 ; triangles qui n''ont qu''on sommet en commun.183 ; Now we have this list, we will make sure that we do not have triangles 184 ; with only a common summit. 189 185 ; 190 186 test = bytarr(nx, ny) … … 192 188 dejavire = 1b-test 193 189 ; 194 tempdeux = systime(1) ; pour key_performance =2190 tempdeux = systime(1) ; For key_performance =2 195 191 vire1 = 0 196 192 vire2 = 0 … … 198 194 vire1 = where( (test*shift(test, -1, -1) $ 199 195 *(1-shift(test, 0, -1))*(1-shift(test, -1, 0))) EQ 1) 200 if vire1[0] NE -1 THEN test[vire1] = 0 ; on vire le rectangle196 if vire1[0] NE -1 THEN test[vire1] = 0 ; We delete the rectangle 201 197 vire2 = where( ((1-test)*(1-shift(test, -1, -1)) $ 202 198 *shift(test, 0, -1)*shift(test, -1, 0)) EQ 1) 203 ; on vire le rectangle du dessus (meme indice x mais egale a1)199 ; We delete the top rectangle (same x index but equal to 1) 204 200 if vire2[0] NE -1 THEN test[vire2+nx] = 0 205 201 ENDWHILE … … 212 208 ; 213 209 if avirer[0] NE -1 then begin 214 tempdeux = systime(1) ; pour key_performance =2210 tempdeux = systime(1) ; For key_performance =2 215 211 indnx = n_elements(listrect) 216 212 indny = n_elements(avirer) -
trunk/SRC/ToBeReviewed/TRIANGULATION/completecointerre.pro
r134 r150 1 ;------------------------------------------------------------2 ;------------------------------------------------------------3 ;------------------------------------------------------------4 ;+5 ; NAME: COMPLETECOINTERRE6 ;7 ; PURPOSE: pour colorier proprement les continents! (c''est une longue8 ; histoire...)9 ;10 ; CATEGORY: pour plt11 ;12 ; CALLING SEQUENCE: completecointerre13 ;14 ; INPUTS: non15 ;16 ; KEYWORD PARAMETERS: _EXTRA17 ;18 ; CONT_COLOR: the color of the continent. defaut value is19 ; (!d.n_colors - 1) < 255 => white20 ;21 ; OUTPUTS: non22 ;23 ; COMMON BLOCKS: common.pro24 ;25 ; SIDE EFFECTS:26 ;27 ; RESTRICTIONS:28 ;29 ; EXAMPLE:30 ;31 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)32 ; 01/10/199933 ;-34 ;------------------------------------------------------------35 ;------------------------------------------------------------36 ;------------------------------------------------------------37 1 PRO draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 38 2 ; … … 55 19 ;------------------------------------------------------------ 56 20 ;------------------------------------------------------------ 21 ;------------------------------------------------------------ 22 ;+ 23 ; 24 ; @file_comments 25 ; To color cleanly continents 26 ; 27 ; @categories 28 ; graphic 29 ; 30 ; @keyword _EXTRA 31 ; Used to pass your keywords 32 ; 33 ; @keyword CONT_COLOR 34 ; The color of the continent. defaut value is 35 ; (!d.n_colors - 1) < 255 => white 36 ; 37 ; @uses 38 ; common.pro 39 ; 40 ; @history 41 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 42 ; 01/10/1999 43 ; 44 ; @version 45 ; $Id$ 46 ; 47 ;- 48 ;------------------------------------------------------------ 49 ;------------------------------------------------------------ 50 ;------------------------------------------------------------ 57 51 PRO completecointerre, COINMONTE = coinmonte, COINDESCEND = coindescend $ 58 52 , CONT_COLOR = cont_color, INDICEZOOM = indicezoom $ … … 66 60 ; if NOT keyword_set(coindescend) then return 67 61 ; if NOT keyword_set(indicezoom) then return 68 tempsun = systime(1) ; pour key_performance62 tempsun = systime(1) ; For key_performance 69 63 ;------------------------------------------------------------ 70 ; definitions des vecteurs coinmont etcoindesc64 ; definitions of vectors coinmont and coindesc 71 65 ;------------------------------------------------------------ 72 66 if keyword_set(coinmonte) then coinmont = coinmonte $ … … 76 70 IF NOT keyword_set(cont_color) THEN cont_color = (!d.n_colors-1) < 255 77 71 ;------------------------------------------------------------ 78 ; definition descoordonnees des points numerotes 1,2,3,4,5,6 cf. les 79 ; schemas en dessous! 72 ; definition of coordinates of points numbered 1,2,3,4,5,6 (see figures below) 80 73 ;------------------------------------------------------------ 81 tempdeux = systime(1) ; pour key_performance =274 tempdeux = systime(1) ; For key_performance =2 82 75 if coinmont[0] NE -1 OR coindesc[0] NE -1 then BEGIN 83 76 if keyword_set(indicezoom) then BEGIN … … 122 115 ; 123 116 ; 124 ; cas coin terre en montee: 125 ; 2 points terre en diagonale montante avec 2 points mer sur 126 ; la diagonale descendante. 117 ; Case land corner in ascent: 118 ; 2 land points in diagonal ascending with 2 ocean points on the descendant diagonal. 127 119 ; 128 120 ; 4 … … 139 131 ; 140 132 if coinmont[0] NE -1 then BEGIN 141 tempdeux = systime(1) ; pour key_performance =2133 tempdeux = systime(1) ; For key_performance =2 142 134 for id = 0, n_elements(coinmont)-1 do BEGIN 143 135 i = coinmont[id] … … 159 151 ENDIF 160 152 ;------------------------------------------------------------ 161 ; cas coin terre en descendante.: 162 ; 2 points terre en diagonale descendante avec 2 points mer sur 163 ; la diagonale montante 164 ; 153 ; Case land corner in descent: 154 ; 2 land points in diagonal descending with 2 ocean points on the ascendant diagonal. 155 ; 165 156 ; 4 166 157 ; t(i+nx)=1 u(i+nx) t(i+nx+1)=0 … … 175 166 ; 176 167 if coindesc[0] NE -1 then begin 177 tempdeux = systime(1) ; pour key_performance =2168 tempdeux = systime(1) ; For key_performance =2 178 169 for id = 0, n_elements(coindesc)-1 do BEGIN 179 170 i = coindesc[id] -
trunk/SRC/ToBeReviewed/TRIANGULATION/definetri.pro
r134 r150 1 1 ;+ 2 ; NAME:definetri3 2 ; 4 ; PURPOSE:Define a triangulation array like TRIANGULATE. 3 ; @file_comments 4 ; efine a triangulation array like TRIANGULATE. 5 5 ; But in a VERY SIMPLE CASE: 6 6 ; the points are regulary-gridded on nx*ny array. … … 35 35 ; 36 36 ; 37 ; CATEGORY: to understand how TRIANGULATE and TRIANGULATION work! 37 ; @categories 38 ; utilities 39 ; 40 ; @param NX {in}{required} 41 ; The x dimension array 38 42 ; 39 ; CALLING SEQUENCE:triangles=definetri(nx, ny [,downward]) 40 ; 41 ; INPUTS: nx and ny are the array dimensions 43 ; @param NY {in}{required} 44 ; The y dimension array 42 45 ; 43 ; OPTIONAL INPUTS: 46 ; @param DOWNWARD {in}{optional} 47 ; When downward is undefine all rectangles are cut 48 ; in using the upward diagonal. Downward is a vector which 49 ; contains the rectangles numbers which are cut in using the 50 ; downward diagonal. 51 ; The rectangle number is define by the index (in a nx*ny 52 ; vector) of the lower-left corner of the rectangle. 44 53 ; 45 ; downward: When downward is undefine all rectangles are cut 46 ; in using the upward diagonal. Downward is a vector which 47 ; contains the rectangles numbers which are cut in using the 48 ; downward diagonal. 49 ; The rectangle number is define by the index (in a nx*ny 50 ; vector) of the lower-left corner of the rectangle. 54 ; @returns 55 ; triangles is a 2d array and is dimensions are 3 and 56 ; 2*(nx-1)*(ny-1) 57 ; triangles is define like in the TRIANGULATE procedure. 51 58 ; 52 ; KEYWORD PARAMETERS: 53 ; 54 ; OUTPUTS: 55 ; triangles is a 2d array and is dimensions are 3 and 56 ; 2*(nx-1)*(ny-1) 57 ; triangles is define like in the TRIANGULATE procedure. 58 ; 59 ; 60 ; OPTIONAL OUTPUTS: 61 ; 62 ; COMMON BLOCKS: 63 ; 64 ; SIDE EFFECTS: 65 ; 66 ; RESTRICTIONS: 67 ; 68 ; PROCEDURE: 69 ; 70 ; EXAMPLE: 59 ; @examples 71 60 ; 72 61 ; triangles=definetri(3,3,[1,3]) … … 84 73 ; 85 74 ; 86 ; MODIFICATION HISTORY: sebastien Masson (smlod@ipsl.jussieu.fr) 75 ; @history 76 ; sebastien Masson (smlod@ipsl.jussieu.fr) 87 77 ; 4/3/1999 88 78 ; 79 ; @version 80 ; $Id$ 89 81 ;- 90 82 FUNCTION definetri, nx, ny, downward -
trunk/SRC/ToBeReviewed/TRIANGULATION/definetri_e.pro
r134 r150 13 13 14 14 ;+ 15 ; NAME:definetri16 15 ; 17 ; PURPOSE:Define a triangulation array like TRIANGULATE but for a 16 ; @file_comments 17 ; Define a triangulation array like TRIANGULATE but for a 18 18 ; E-grid type 19 19 ; 20 ; CATEGORY: make contours with E-grid type 20 ; @categories 21 ; Make contours with E-grid type 22 ; 23 ; @param NX {in}{required} 24 ; The x dimension array 21 25 ; 22 ; CALLING SEQUENCE:triangles=definetri(nx, ny [,vertical]) 23 ; 24 ; INPUTS: nx and ny are the array dimensions 26 ; @param NY {in}{required} 27 ; The y dimension array 25 28 ; 26 ; OPTIONAL INPUTS: 29 ; @param SINGULAR {in}{optional} 30 ; When singular is undefined all rectangles are cut 31 ; in using the vertical diagonal. Singular is a vector which 32 ; contains the rectangles numbers which are cut in using the 33 ; horizontal diagonal. 34 ; The rectangle number is define by the index (in a nx*ny 35 ; vector) of the lower-left corner of the rectangle. 27 36 ; 28 ; vertical: When vertical is undefine all rectangles are cut 29 ; in using the horizontal diagonal. Vertical is a vector which 30 ; contains the rectangles numbers which are cut in using the 31 ; vertical diagonal. 32 ; The rectangle number is define by the index (in a nx*ny 33 ; vector) of the lower-left corner of the rectangle. 34 ; 35 ; KEYWORD PARAMETERS: 36 ; 37 ; OUTPUTS: 38 ; triangles is a 2d array and is dimensions are 3 and 39 ; 2*(nx-1)*(ny-1) 40 ; triangles is define like in the TRIANGULATE procedure. 41 ; 42 ; 43 ; OPTIONAL OUTPUTS: 44 ; 45 ; COMMON BLOCKS: 46 ; 47 ; SIDE EFFECTS: 48 ; 49 ; RESTRICTIONS: 50 ; 51 ; PROCEDURE: 52 ; 53 ; EXAMPLE: 37 ; @keyword SHIFTED 54 38 ; 55 39 ; 56 ; MODIFICATION HISTORY: sebastien Masson (smlod@ipsl.jussieu.fr) 40 ; @returns 41 ; Triangles is a 2d array and is dimensions are 3 and 42 ; (nx-1)*(ny-1) 43 ; Triangles is defined like in the TRIANGULATE procedure. 44 ; 45 ; @history 46 ; Sebastien Masson (smlod@ipsl.jussieu.fr) 57 47 ; June 2001 58 48 ; 49 ; @version 50 ; $Id$ 51 ; 52 ; @todo 53 ; seb: documenter SHIFTED 59 54 ;- 60 55 FUNCTION definetri_e, nx, ny, singular, SHIFTED = shifted -
trunk/SRC/ToBeReviewed/TRIANGULATION/dessinetri.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME:dessinetri6 5 ; 7 ; PURPOSE:dessine la triangulation 6 ; @file_comments 7 ; Draw the triangulation 8 8 ; 9 ; CATEGORY:pour comprendre comment ca marche 9 ; @categories 10 ; utilities 11 ; 12 ; @param TRI {in}{optional} 13 ; Array which define the triangulation (provided by triangule.pro or triangulate) 14 ; 15 ; @param X {in}{optional} 16 ; The x position of points to which the trangulation refer to 17 ; (see the x array provided in triangulate) 10 18 ; 11 ; CALLING SEQUENCE:dessinetri [, tri, x, y] 12 ; 13 ; INPUTS:optionnels 14 ; par defaut on choisit la triangulation qui est utilise pour 15 ; les plots et on la trace aux points definites par vargrid 19 ; @param Y {in}{optional} 20 ; The y position of points to which the trangulation refer to 21 ; (see the y array provided in triangulate) 16 22 ; 17 ; sinon il faut fournir les tableaux 18 ; tri definissant la triangulation (fournis par triangule.pro 19 ; ou triangulate) 20 ; x et y qui sont les positions de points a laquelle se raporte 21 ; la triangulation (cf. les tableau x et y fournis ds 22 ; triangulate) 23 ; @keyword WAIT 24 ; =x. to call wait x second between each triangle draw. 23 25 ; 24 ; KEYWORD PARAMETERS: 26 ; @keyword ONEBYONE 27 ; To draw the triangles one by one 25 28 ; 26 ; All plots or polyfill keywords. 29 ; @keyword FILL 30 ; To fill the triangles (using polyfill) instead of plotting them 27 31 ; 28 ; WAIT=x. to call wait x second between each triangle draw. 32 ; @keyword CHANGECOLOR 33 ; =n. To change the color of each traingle. n colors 34 ; will be used and repeted if necessary. 29 35 ; 30 ; /ONEBYONE: to draw the triangles one by one 36 ; @uses 37 ; common.pro 31 38 ; 32 ; /FILL: to fill the triangles (using polyfill) instead of plotting them 39 ; @history 40 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 33 41 ; 34 ; CHANGECOLOR=n. to change the color of each traingle. n colors 35 ; will be used and repeted if necessary. 36 ; 37 ; 38 ; OUTPUTS: 39 ; 40 ; COMMON BLOCKS:common.pro 41 ; 42 ; SIDE EFFECTS: 43 ; 44 ; RESTRICTIONS: 45 ; 46 ; EXAMPLE: 47 ; 48 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) 42 ; @version 43 ; $Id$ 49 44 ; 50 45 ;- … … 58 53 ; 59 54 @common 60 tempsun = systime(1) ; pour key_performance55 tempsun = systime(1) ; For key_performance 61 56 a = '' 62 57 if n_params() EQ 3 then BEGIN … … 99 94 color = color#replicate(1, n_elements(tri)/3/n_elements(color)+1) 100 95 ; 101 tempdeux = systime(1) ; pour key_performance =296 tempdeux = systime(1) ; For key_performance =2 102 97 for i = 0L, n_elements(tri)/3-1 do begin 103 98 t = [tri[*, i], tri[0, i]] -
trunk/SRC/ToBeReviewed/TRIANGULATION/drawcoast_c.pro
r134 r150 1 ;+ 2 ; @file_comments 3 ; 4 ; 5 ; @categories 6 ; 7 ; 8 ; @param MASK 9 ; 10 ; 11 ; @param XF 12 ; 13 ; 14 ; @param YF 15 ; 16 ; 17 ; @param NX 18 ; 19 ; 20 ; @param NY 21 ; 22 ; 23 ; @keyword COAST_COLOR 24 ; The color of the coastline. 25 ; defaut value is 0 => black 26 ; 27 ; @keyword COAST_THICK 28 ; The thick of the trait to trace continents 29 ; By default, it is 1. 30 ; 31 ; @keyword XSEUIL 32 ; To eliminate segments of coasts which are to big (which link points which can 33 ; be close on the sphere but distant on the drawing). We delete all segments 34 ; whose the size surpass the size of the window following X/XSEUIL. By default, 35 ; XSEUIL=5. But it can be to big if we do a big zoom or a little one for some 36 ; projections... We specify it the keyword thanks to this keyword! 37 ; 38 ; @keyword YSEUIL 39 ; See XSEUIL 40 ; 41 ; @keyword _EXTRA 42 ; Used to pass our keywords 43 ; 44 ; @returns 45 ; 46 ; 47 ; @uses 48 ; 49 ; 50 ; @restrictions 51 ; 52 ; 53 ; @examples 54 ; 55 ; 56 ; @history 57 ; 58 ; 59 ; @version 60 ; 61 ; @todo 62 ; Seb: remplir le header 63 ; 64 ;- 65 ;--------------------------------------------------------- 1 66 PRO drawcoast_c, mask, xf, yf, nx, ny, COAST_COLOR = coast_color, COAST_THICK = coast_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 2 ;---------------------------------------------------------3 67 ; 4 68 compile_opt idl2, strictarrsubs … … 10 74 ENDIF 11 75 ;--------------------------------------------------------- 12 tempsun = systime(1) ; pour key_performance76 tempsun = systime(1) ; For key_performance 13 77 ;--------------------------------------------------------- 14 78 ; 15 ; on trace les segments verticaux:79 ; We trace vertical segments: 16 80 ; 17 81 if NOT keyword_set(yseuil) then yseuil = 5. < (min([nx, ny])-2) 18 82 distanceseuil = (!p.position[3]-!p.position[1])/yseuil 19 ; liste: liste des points i pourlesquels on va tracer un segment entre 20 ; le point i,j-1 et i,j 21 tempdeux = systime(1) ; pour key_performance =2 83 ; list: list of points i for which we will trace a segment between the point i,j-1 and i,j 84 tempdeux = systime(1) ; For key_performance =2 22 85 liste = where((mask+shift(mask, -1, 0)) EQ 1 $ 23 86 AND ((xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2) LE distanceseuil^2) 24 87 IF liste[0] NE -1 THEN BEGIN 25 ; on recupere lx et ly qui sont les indices ds un tableau 2d des 26 ; points donnes par liste 88 ; We recuperate lx an dly which are indexes in a 2d array of points given by list 27 89 ly = liste/nx & lx = temporary(liste)-nx*ly 28 indice = where(ly NE 0) ; on ne prend pas les points concernant90 indice = where(ly NE 0) ; We do not take points concerning 29 91 if indice[0] NE -1 then begin 30 ; la premiere ligne car ds ce cas le pt j-1 n''est pas definit92 ; the first line because in this case, the point j-1 is undefined. 31 93 lx = lx[indice] & ly = ly[temporary(indice)] 32 ; boucle sur les points concernes et trace du segment33 ; rq: on utilise plost au lieu de plot car plots est bcp plus rapide.94 ; Loop on concerned points and drawing of the segment. 95 ; Comment: we use plots instead of plot because plots goes faster. 34 96 IF testvar(var = key_performance) EQ 2 THEN $ 35 97 print, 'temps tracecote: determiner liste des points concernes par un trait vertical', systime(1)-tempdeux 36 tempdeux = systime(1) ; pour key_performance =298 tempdeux = systime(1) ; For key_performance =2 37 99 for pt = 0L, n_elements(lx)-1 do BEGIN 38 100 i = lx[pt] & j = ly[pt] … … 45 107 ENDIF 46 108 ; 47 ; pour le trace des segments horizontaux, c''est la meme chose sauf 48 ; qu'il faut faire attention si on est periodique: 109 ; For the drawing of horizontal segments , it is the same thing but we have to be careful if it is periodic. 49 110 ; 50 ; si on est periodique on duplique la premiere colonne et on la met a 51 ; la fin. (ceci est fait non pas pour le shift qui est par defaut 52 ; periodique mais pour le plots 53 tempdeux = systime(1) ; pour key_performance =2 111 ; If it is periodic, we duplicate the first column and we put it at the end. 112 ; (This is made not for the shift which is periodic by default but for the plots) 113 tempdeux = systime(1) ; For key_performance =2 54 114 if keyword_set(key_periodic) AND nx EQ jpi then begin 55 115 mask = [mask, mask[0, *]] … … 66 126 indice = where(ly NE ny-1 AND lx NE 0) 67 127 if indice[0] NE -1 then begin 68 ; on ne prend pas les points de la 69 ; premiere colonne et de la derniere ligne (car on l''a rajoute artificiellement!)) 128 ; We do not take points of the first column and of the last line (because we added artificially) 70 129 lx = lx[indice] & ly = ly[temporary(indice)] 71 130 IF testvar(var = key_performance) EQ 2 THEN $ 72 131 print, 'temps tracecote: determiner liste des points concernes par un trait horizontal', systime(1)-tempdeux 73 tempdeux = systime(1) ; pour key_performance =2132 tempdeux = systime(1) ; For key_performance =2 74 133 for pt = 0L, n_elements(lx)-1 do BEGIN 75 134 i = lx[pt] & j = ly[pt] -
trunk/SRC/ToBeReviewed/TRIANGULATION/drawcoast_e.pro
r134 r150 1 ;+ 2 ; @file_comments 3 ; 4 ; 5 ; @categories 6 ; 7 ; 8 ; @param MASK 9 ; 10 ; 11 ; @param XF 12 ; 13 ; 14 ; @param YF 15 ; 16 ; 17 ; @param NX 18 ; 19 ; 20 ; @param NY 21 ; 22 ; 23 ; @keyword COAST_COLOR 24 ; The color of the coastline. 25 ; defaut value is 0 => black 26 ; 27 ; @keyword COAST_THICK 28 ; The thick of the trait to trace continents 29 ; By default, it is 1. 30 ; 31 ; @keyword XSEUIL 32 ; To eliminate segments of coasts which are to big (which link points which can 33 ; be close on the sphere but distant on the drawing). We delete all segments 34 ; whose the size surpass the size of the window following X/XSEUIL. By default, 35 ; XSEUIL=5. But it can be to big if we do a big zoom or a little one for some 36 ; projections... We specify it the keyword thanks to this keyword! 37 ; 38 ; @keyword YSEUIL 39 ; See XSEUIL 40 ; 41 ; @keyword _EXTRA 42 ; Used to pass our keywords 43 ; 44 ; @returns 45 ; 46 ; 47 ; @uses 48 ; 49 ; 50 ; @restrictions 51 ; 52 ; 53 ; @examples 54 ; 55 ; 56 ; @history 57 ; 58 ; 59 ; @version 60 ; 61 ; @todo 62 ; Seb: remplir le header 63 ; 64 ;- 65 1 66 PRO drawcoast_e, mask, xf, yf, nx, ny, COAST_COLOR = coast_color, COAST_THICK = coast_thick, YSEUIL = yseuil, XSEUIL = xseuil, onemore = onemore, _extra = ex 2 67 ;--------------------------------------------------------- -
trunk/SRC/ToBeReviewed/TRIANGULATION/drawsectionbottom.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME:drawsectionbottom6 5 ; 7 ; PURPOSE:fill and draw the bottom continents for a real section. 6 ; @file_comments 7 ; Fill and draw the bottom continents for a real section. 8 8 ; 9 ; CATEGORY: 9 ; @categories 10 ; 11 ; @param MASKIN {in}{requierd} 10 12 ; 11 ; CALLING SEQUENCE: 12 ; 13 ; INPUTS: 13 ; @param XXAXISIN {in}{requierd} 14 14 ; 15 ; KEYWORD PARAMETERS:15 ; @param DEPTHSIN {in}{requierd} 16 16 ; 17 ; COAST_COLOR: the color of the coastline. 18 ; defaut value is 0 => black 17 ; @keyword COAST_COLOR 18 ; The color of the coastline. 19 ; defaut value is 0 => black 19 20 ; 20 ; COAST_THICK: the thickness of the coastline. 21 ; defaut value is 1 21 ; @keyword COAST_THICK 22 ; The thickness of the coastline. 23 ; defaut value is 1 22 24 ; 23 ; CONT_COLOR: the color of the continent. defaut value is 24 ; (!d.n_colors - 1) < 255 => white 25 ; @keyword CONT_COLOR 26 ; The color of the continent. defaut value is 27 ; (!d.n_colors - 1) < 255 => white 25 28 ; 26 ; OUTPUTS: 29 ; @uses 30 ; common.pro 27 31 ; 28 ; COMMON BLOCKS:common.pro 29 ; 30 ; SIDE EFFECTS: 31 ; 32 ; RESTRICTIONS:simple way to fill continents for a section (using the 32 ; @restrictions 33 ; Simple way to fill continents for a section (using the 33 34 ; fact that continents are wider at the bottom than at the top). 34 35 ; 35 ; EXAMPLE: 36 ; @history 37 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 38 ; June 14, 2002 36 39 ; 37 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) 38 ; June 14, 2002 40 ; @version 41 ; $Id$ 42 ; 43 ; @todo 44 ; Seb: definir params 39 45 ;- 40 46 ;------------------------------------------------------------ -
trunk/SRC/ToBeReviewed/TRIANGULATION/fillcornermask.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME: FILLCORNERMASK6 5 ; 7 ; PURPOSE: pour colorier proprement les continents! (c''est une longue8 ; histoire...)6 ; @file_comments 7 ; To color cleanly continents 9 8 ; 10 ; CATEGORY: pour plt 9 ; @categories 10 ; graphic 11 11 ; 12 ; CALLING SEQUENCE: completecointerre 13 ; 14 ; INPUTS: non 12 ; @keyword _EXTRA 13 ; Used to pass your keywords 15 14 ; 16 ; KEYWORD PARAMETERS: _EXTRA 15 ; @keyword CONT_COLOR 16 ; The color of the continent. defaut value is 17 ; (!d.n_colors - 1) < 255 => white 17 18 ; 18 ; CONT_COLOR: the color of the continent. defaut value is19 ; (!d.n_colors - 1) < 255 => white19 ; @uses 20 ; common.pro 20 21 ; 21 ; OUTPUTS: non 22 ; @history 23 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 24 ; 8/8/2002 22 25 ; 23 ; COMMON BLOCKS: common.pro 26 ; @version 27 ; $Id$ 24 28 ; 25 ; SIDE EFFECTS:26 ;27 ; RESTRICTIONS:28 ;29 ; EXAMPLE:30 ;31 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)32 ; 8/8/200233 29 ;- 34 30 ;------------------------------------------------------------ … … 45 41 if NOT keyword_set(coinmonte) AND NOT keyword_set(coindescend) then return 46 42 ; 47 tempsun = systime(1) ; pour key_performance43 tempsun = systime(1) ; For key_performance 48 44 ; 49 45 IF NOT keyword_set(cont_color) THEN cont_color = (!d.n_colors-1) < 255 50 46 ;------------------------------------------------------------ 51 ; definition descoordonnees des points numerotes 1,2,3,4,5,6 cf. les 52 ; schemas en dessous! 47 ; definition of coordinates of points numbered 1,2,3,4,5,6 (see figures below) 53 48 ;------------------------------------------------------------ 54 49 ; … … 64 59 ; 65 60 ; 66 ; cas coin terre en montee: 67 ; 2 points terre en diagonale montante avec 2 points mer sur 68 ; la diagonale descendante. 61 ; Case land corner in ascent: 62 ; 2 land points in diagonal ascending with 2 ocean points on the descendant diagonal. 69 63 ; 70 64 ; 3 … … 96 90 endif 97 91 ;------------------------------------------------------------ 98 ; cas coin terre en descendante.: 99 ; 2 points terre en diagonale descendante avec 2 points mer sur 100 ; la diagonale montante 92 ; Case land corner in descent: 93 ; 2 land points in diagonal descending with 2 ocean points on the ascendant diagonal. 101 94 ; 102 95 ; 4 -
trunk/SRC/ToBeReviewed/TRIANGULATION/section.pro
r134 r150 1 ;------------------------------------------------------------2 ;------------------------------------------------------------3 ;------------------------------------------------------------4 1 ;+ 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 ; 2 ; @file_comments 3 ; 4 ; 5 ; @categories 6 ; 7 ; 8 ; @param FIELD 9 ; 10 ; 11 ; @param RES 12 ; 13 ; 14 ; @param GLAMAXE 15 ; 16 ; 17 ; @param GPHIAXE 18 ; 19 ; 20 ; @keyword ENDPOINTS 21 ; 22 ; 23 ; @keyword BOXZOOM 24 ; 25 ; 26 ; @keyword TYPE 27 ; 28 ; 29 ; @keyword WDEPTH 30 ; 31 ; 32 ; @keyword DIREC 33 ; 34 ; 35 ; @keyword SHOWBUILD 36 ; 37 ; 38 ; @keyword ONLYBOX 39 ; 40 ; 41 ; @keyword _EXTRA 42 ; Used to pass your keywords 43 ; 44 ; @returns 45 ; 46 ; 47 ; @uses 48 ; common.pro 49 ; 50 ; @restrictions 51 ; 52 ; 53 ; @examples 54 ; 55 ; 56 ; @history 57 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 58 ; 59 ; @version 60 ; 29 61 ;- 30 ;------------------------------------------------------------31 ;------------------------------------------------------------32 ;------------------------------------------------------------33 62 34 63 PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints $ … … 52 81 ;--------------------- 53 82 ;--------------------- 54 ; definition de boxzoom en fonction de endpoints 55 ; puis redefinition du domaine 83 ; definition of boxzoom in function of endpoints, then redefinition of the domain 56 84 boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $ 57 85 , min([endpoints[1], endpoints[3]], max = ma13), ma13] … … 83 111 ; how can we find a good test value? 84 112 IF nx * ny LE 10 THEN domdef, GRIDTYPE = grillechoice, _extra = ex 85 ; on redefinit lon1, ... au cas ou findalways ait ete utilise dsdomdef113 ; We redefine lon1, ... in case findalways has been used in domdef 86 114 lon1 = min([endpoints[0], endpoints[2]], max = lon2) 87 115 lat1 = min([endpoints[1], endpoints[3]], max = lat2) … … 89 117 ; until its bottom part. 90 118 if strpos(type, 'z') NE -1 THEN BEGIN 91 ; on garde les yranges (axe z) avant de changer laboxzoom.119 ; We keep yranges (axis z) before changing the boxzoom. 92 120 !y.range = [localbox[nelbox-1], localbox[nelbox-2]] 93 121 if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN … … 105 133 ENDIF 106 134 !y.range = [localbox[nelbox-1], localbox[nelbox-2]] 107 ; on recupere la grille sur la boxzoom135 ; We recuperate the grid on the boxzoom. 108 136 ; 109 137 IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat' … … 126 154 127 155 ;--------------------- 128 ; on definit la triangulation qui va nous permetre de determiner la129 ; section. on la recalcule car elle doit etre definie sur la terre130 ; aussi bien que sur la mer. suivant le sens de la section -plutot131 ; longitude ou plutot latitude- on definit la facon de trianguler156 ; We define the triangulation which will allows us to determinate the section. 157 ; We recalculate it because it must be defined on the Earth and on oceans. 158 ; Following the direction of the section -rather longitude or rather latitude- 159 ; we define the way to triangulate. 132 160 if strpos(type, 'x') NE -1 then BEGIN 133 161 downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] … … 146 174 ENDIF 147 175 ;--------------------- 148 ; equation de la droite suivant laquelle on fait la section176 ; Equation of the line on which we do the section. 149 177 ;--------------------- 150 178 abc = linearequation(endpoints[0:1], endpoints[2:3]) 151 179 glamtri = glam[tri] 152 180 gphitri = gphi[tri] 153 ; quels sont les points de la triangulation qui sont au dessus et au 154 ; dessous de la droite ? 181 ; Which points of the triangulation are above and below the line? 155 182 if abc[1] NE 0 THEN $ 156 183 test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $ … … 158 185 159 186 zero123 = total(test, 1) 160 ; to keep: triangles de la triangulation qui sont a cheval sur la droite.187 ; to keep: triangles of the triangulation which are over the line. 161 188 tokeep1 = where(zero123 EQ 1) 162 189 tokeep2 = where(temporary(zero123) EQ 2) … … 165 192 test = test[*, tokeep] 166 193 tri = tri[*, tokeep] 167 ; quel est le sommet du triangle qui est seul d''un cote de la droite?194 ; Which summit of the triangle is alone in a side of the line? 168 195 single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1) 169 196 single1 = single1-(single1/3)*3 … … 177 204 178 205 single = [temporary(single1), temporary(single2)] 179 ; points1 le point du triangles qui est seul d''un cote de la droite.180 ; point2 l''autre point du triangle de l''autre cote de la droite206 ; points1 the point (of the triangle) alone in a side of the line. 207 ; point2 the other point of the triangle in the other side of the line. 181 208 point1 = [single, single] 182 209 point2 = [single EQ 0, 1 + (single LE 1)] … … 189 216 points1 = tri[point1, index] 190 217 points2 = tri[point2, temporary(index)] 191 ; points : complex e contenant les couples de points de part et192 ; d''autre de la droite. Ils faut supprimer les doublons.218 ; points : complex containing couples of points in a side and the other 219 ; side of the line. We have to delete duplicates. 193 220 points = dcomplex(points1, points2) 194 221 points = points[uniq(points, sort(points))] 195 222 symetrique = dcomplex(imaginary(points), double(points)) 196 223 points = points[where(points-shift(temporary(symetrique), 1) NE 0)] 197 ; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite.198 ; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite224 ; points1 coordinates of the point of the triangle which is alone in a side of the line. 225 ; point2 coordinates of the other point of the triangle in the other side of the line. 199 226 points1 = complex(glam[ double(points)], gphi[ double(points)]) 200 227 points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) 201 ; droites les equations des droites dont on cherche l''intersection 202 ; avec la section. 228 ; droites equations of line whose we look for the intersection wit the section. 203 229 droites = linearequation(points1, points2) 204 230 inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) 205 231 206 ; les ccordonnes geographiques des points que l''on cherche sur lasection.232 ; Geographic coordinates of points we look for on the section. 207 233 glamaxe = float(inter) 208 234 gphiaxe = imaginary(inter) 209 ; on les range ds l''ordre croissant entre les bornes de la section235 ; We arrange them in the growing order between boundaries of the section. 210 236 if strpos(type, 'x') NE -1 then BEGIN 211 237 sort = sort(glamaxe) … … 272 298 poids = poids#replicate(1, nz) 273 299 res = poids*value1+(1-poids)*value2 274 ; moyenne suivantz ?300 ; average following z ? 275 301 if strpos(type, 'z') EQ -1 then begin 276 302 nan = where(finite(res) EQ 0) … … 319 345 poids = reform(poids, npoints, nz, jpt, /over) 320 346 res = poids*value1+(1-poids)*value2 321 ; moyenne suivantz ?347 ; average following z ? 322 348 if strpos(type, 'z') EQ -1 then begin 323 349 nan = where(finite(res) EQ 0) -
trunk/SRC/ToBeReviewed/TRIANGULATION/tracecote.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME:tracecote6 5 ; 7 ; PURPOSE: dessine les cotes ds plt 6 ; @file_comments 7 ; Draw coasts in plt. 8 8 ; 9 ; CATEGORY: pour faire un joli dessin 9 ; @categories 10 ; graphic 10 11 ; 11 ; CALLING SEQUENCE:tracecote,mask12 12 ; 13 13 ; INPUTS:mask le tableau mask sur la zone consideree pour le dessin 14 14 ; 15 ; KEYWORD PARAMETERS: 15 ; @keyword SURFACE_COASTLINE 16 ; To draw the surface coast line instead of 17 ; the coast line at level firstz[tw]. Usefull only for deep 18 ; plots! 16 19 ; 17 ; COAST_COLOR: the color of the coastline.18 ; defaut value is 0 => black20 ; @keyword _EXTRA 21 ; used to pass your keywords 19 22 ; 20 ; COAST_THICK: l''epaisseur du trait pour tracer les21 ; continents. par defaut c''est 1.23 ; @uses 24 ; common.pro 22 25 ; 23 ; /SURFACE_COASTLINE: to draw the furface coast line instead of24 ; the coast line at level firstz[tw]. Usefull only for deep25 ; plots!26 ; @history 27 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 28 ; 30/9/1999 26 29 ; 27 ; XSEUIL: pour eliminer les segments de cote qui sont trop 28 ; grand (qui relient des points qui peuvent etre tres proches 29 ; sur la sphere mais tres eloignes sur le dessin) on supprime 30 ; tous les egments dot la taille depasse: 31 ; taille de la fenetre suivant X/ xseuil. 32 ; Par defaut xseuil est egale a 5. masi peut etre trop grand si 33 ; on fait un fort zoom ou trout petit pour certaines 34 ; projections... le specifier alors a l''aide de ce mot cle! 30 ; @version 31 ; $Id$ 35 32 ; 36 ; YSEUIL: cf. xseuil37 ;38 ; OUTPUTS: rien39 ;40 ; COMMON BLOCKS:common.pro41 ;42 ; SIDE EFFECTS:43 ;44 ; RESTRICTIONS:45 ;46 ; EXAMPLE:47 ;48 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)49 ; 30/9/199950 33 ;- 51 34 ;------------------------------------------------------------ … … 64 47 ENDIF 65 48 ;-------------------------------------------------------------- 66 tempsun = systime(1) ; pour key_performance49 tempsun = systime(1) ; For key_performance 67 50 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 68 51 ; 69 ; on agrandi un peu le cadre definit par les premier..., dernier... de 70 ; facon a bien recuperer les bords de cote qui sont en bordure du 71 ; domaine a tracer 52 ; We enlarge a bit the frame defined by firsts..., lasts... in order to 53 ; recuperate edges of the coast which are in the edging of the domain. 72 54 ; 73 tempdeux = systime(1) ; pour key_performance =255 tempdeux = systime(1) ; For key_performance =2 74 56 firstx = 0 > (min([firstxt, firstxf])-1) 75 57 lastx = (max([lastxt, lastxf])+1) < (jpi-1) … … 78 60 nx = lastx-firstx+1 79 61 ny = lasty-firsty+1 80 ; quel niveau vertical choisir?62 ; Which vertical level choose? 81 63 IF keyword_set(surface_coastline) THEN firstz = 0 ELSE $ 82 64 IF strupcase(vargrid) eq 'W' THEN firstz = firstzw ELSE firstz = firstzt 83 ; attribution du masque et des coordonnes delimitant les limites de la 84 ; terre (coordonnees f) 65 ; Attribution of the mask and of coordinates delimiting limits of the land (coordinates f). 85 66 mask = tmask[firstx:lastx, firsty:lasty, firstz] 86 67 xf = glamf[firstx:lastx, firsty:lasty] … … 90 71 ; 91 72 if key_gridtype EQ 'e' then onemore = xf[0, 0] gT xf[0, 1] 92 ; on passe en coordonnee normaliser pour pouvoir s'affranchir du type 93 ; de projection choisie et du suport surlequel on fait le dessin 94 ; (ecran ou postscript) 73 ; We pass in normalized coordinates to be able to become independant from the projection's 74 ; type choosen and from the support on which we do the drawing (screen or postscript) 95 75 z = convert_coord(xf[*],yf[*],/data,/to_normal) 96 76 xf = reform(z[0, *], nx, ny) … … 98 78 tempvar = SIZE(TEMPORARY(z)) 99 79 ; 100 ; attention, suivant la projection certains points x ou y peuvent 101 ; devenir NaN (cf. points deriere la terre ds une projection 102 ; orthographique) 80 ; Beware, following the projection, some points x or y can become NaN (see point 81 ; behind the earth in an orthographic projection). 103 82 ; 104 ; on met les points a eliminer a une tres gande valeur comme ca il ne105 ; passerons pas le test avec distanceseuil (cf. plus bas)83 ; We put points to be eliminated at a very big value so that they will not pass the 84 ; test with distanceseuil (see further). 106 85 ; 107 86 if (!map.projection LE 7 AND !map.projection NE 0) $ … … 117 96 ind = where(yf LT !p.position[1] OR yf GT !p.position[3]) 118 97 IF ind[0] NE -1 THEN yf[ind] = 1e5 119 tempvar = SIZE(TEMPORARY(ind)) ; on efface ind98 tempvar = SIZE(TEMPORARY(ind)) ; we delete ind 120 99 ; 121 100 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' -
trunk/SRC/ToBeReviewed/TRIANGULATION/tracemask.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME:tracemask6 5 ; 7 ; PURPOSE:dessiner des contour d''un mask 6 ; @file_comments 7 ; Draw contours of a mask 8 8 ; 9 ; CATEGORY:plus simple que tracecote, car ne s''occuppe pas du type de 10 ; projection et de la periodicite de la grille 9 ; @categories 10 ; utilities 11 ; 12 ; @param MASKENTREE {in}{required} 13 ; 2d array specifing the mask 14 ; 15 ; @param XIN {in}{required}, 16 ; 2d array specifing longitude coordinates. 17 ; 18 ; @param YIN {in}{required}, 19 ; 2d array specifing latitude coordinates. 11 20 ; 12 ; CALLING SEQUENCE: tracemask, maskentree, xentree, yentree 13 ; 14 ; INPUTS:maskentree, xentree, yentree tableaux 2d specifiant le mask 15 ; et ses coordonees en longitude te latitude. 21 ; @keyword COAST_COLOR 22 ; The color of the coastline. 23 ; defaut value is 0 => black 16 24 ; 17 ; KEYWORD PARAMETERS: 25 ; @keyword COAST_THICK 26 ; The thick of the trait to trace continents 27 ; It is 1 by default. 18 28 ; 19 ; COAST_COLOR: the color of the coastline.20 ; defaut value is 0 => black29 ; @keyword OVERPLOT 30 ; To do a plot over an other one. 21 31 ; 22 ; COAST_THICK: l''epaisseur du trait pour tracer les23 ; continents. par defaut c''est 1.32 ; @keyword _EXTRA 33 ; used to pass your keywords 24 34 ; 25 ; OUTPUTS: none 35 ; @uses 36 ; common.pro 26 37 ; 27 ; COMMON BLOCKS:common.pro 38 ; @history 39 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 28 40 ; 29 ; SIDE EFFECTS: 30 ; 31 ; RESTRICTIONS: 32 ; 33 ; EXAMPLE: 34 ; 35 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) 41 ; @version 42 ; $Id$ 36 43 ; 37 44 ;- … … 51 58 ENDIF 52 59 ;--------------------------------------------------------- 53 tempsun = systime(1) ; pour key_performance54 ; on s''afranchit des problemes de bord:55 tempdeux = systime(1) ; pour key_performance =260 tempsun = systime(1) ; For key_performance 61 ; We avoid edging problems: 62 tempdeux = systime(1) ; For key_performance =2 56 63 tailleentree = size(maskentree) 57 64 nx = tailleentree[1]+1 … … 62 69 IF n_elements(yin) EQ 0 THEN yentree = findgen(ny-1) ELSE yentree = yin 63 70 IF (size(yentree))[0] EQ 1 THEN yentree = replicate(1,nx-1)#yentree 64 ; on agrandi le mask de une colonne a gauche et de une colonne en bas71 ; We enlarge the mask by 1 column to the left an d1 line to the bottom 65 72 mask = intarr(tailleentree[1]+1, tailleentree[2]+1) 66 73 mask[1:tailleentree[1], 1:tailleentree[2]] = maskentree 67 ; les 2 premieres colonnes sont identiques74 ; The 2 first columns are identical. 68 75 mask[0, 1:tailleentree[2]] = maskentree[0, *] 69 ; les 2 premieres lignes sont identiques76 ; The 2 first lines are identical. 70 77 mask[1:tailleentree[1], 0] = maskentree[*, 0] 71 ; on calcul la position suivant x des points qui seviront a tracer le 72 ; masque. ils sont situes entre chaque points du masque, sauf pour la 73 ; derniere colonne que l''on ne peut pas calculer et que l''on met 74 ; donc a max(!x.range) 75 xrange = !x.range[sort(!x.range)] ; si reverse_x est utilise! 78 ; We calculate the position following x of points which will serve to trace the mask. They are situated between each points of the mask, exept for the last column we can not calculate and so we put at max (!x.range). 79 xrange = !x.range[sort(!x.range)] ; if REVERSE_X is used 76 80 xentree = .5*(xentree+shift(xentree, -1, 0)) 77 81 IF not keyword_set(overplot) THEN xentree[nx-2, *] = xrange[1] $ 78 82 ELSE xentree[nx-2, *] = xentree[nx-3, *] 79 ; on seuil83 ; we sill 80 84 xentree = xrange[0] > xentree < xrange[1] 81 ; on agrandit le tableau85 ; we enlarge the array 82 86 xf = fltarr(nx, ny) 83 87 xf[1:nx-1, 1:ny-1] = xentree … … 103 107 print, 'temps tracemask: determination du mask et des ses coordonnes', systime(1)-tempdeux 104 108 ; 105 ; on trace les segments verticaux:109 ; We trace vertical segments: 106 110 ; 107 tempdeux = systime(1) ; pour key_performance =2111 tempdeux = systime(1) ; For key_performance =2 108 112 liste = where(mask+shift(mask, -1, 0) EQ 1) 109 113 IF liste[0] NE -1 THEN BEGIN 110 ; on recupere lx et ly qui sont les indices ds un tableau 2d des 111 ; points donnes par liste 114 ; We recuperate lx and ly which are indexes in a 2d array of points given by list 112 115 ly = liste/nx & lx = temporary(liste)-nx*ly 113 indice = where(ly NE 0) ; on ne prend pas les points concernant114 ; la premiere ligne car ds ce cas le pt j-1 n''est pas definit116 indice = where(ly NE 0) ; We do not take points concernining 117 ; the first line because in this case, the point j-1 is not defined 115 118 if indice[0] NE -1 then begin 116 119 lx = lx[indice] & ly = ly[temporary(indice)] 117 120 IF testvar(var = key_performance) EQ 2 THEN $ 118 121 print, 'temps tracemask: liste traits verticaux', systime(1)-tempdeux 119 tempdeux = systime(1) ; pour key_performance =2120 ; boucle sur les points concernes et trace du segment121 ; rq: on utilise plots au lieu de plot car plots est bcp plus rapide.122 tempdeux = systime(1) ; For key_performance =2 123 ; loop on concerned points and drawing of the segment. 124 ; comments: we use plots instead of plot because plots is faster. 122 125 for pt = 0L, n_elements(lx)-1 do BEGIN 123 126 i = lx[pt] & j = ly[pt] … … 132 135 ENDIF 133 136 ; 134 ; on trace les segments horizontaux:137 ; We trace horizontal segments: 135 138 ; 136 tempdeux = systime(1) ; pour key_performance =2139 tempdeux = systime(1) ; For key_performance =2 137 140 liste = where(mask+shift(mask, 0, -1) EQ 1) 138 141 IF liste[0] NE -1 THEN BEGIN 139 142 ly = liste/nx & lx = temporary(liste)-nx*ly 140 indice = where(lx NE 0) ; on ne prend pas les points de la premiere colonne143 indice = where(lx NE 0) ; We do not take point sof the first column. 141 144 if indice[0] EQ -1 then return 142 145 lx = lx[indice] & ly = ly[temporary(indice)] 143 146 IF testvar(var = key_performance) EQ 2 THEN $ 144 147 print, 'temps tracemask: liste traits horizontaux', systime(1)-tempdeux 145 tempdeux = systime(1) ; pour key_performance =2148 tempdeux = systime(1) ; For key_performance =2 146 149 for pt = 0L, n_elements(lx)-1 do BEGIN 147 150 i = lx[pt] & j = ly[pt] -
trunk/SRC/ToBeReviewed/TRIANGULATION/triangule.pro
r134 r150 1 ;------------------------------------ 2 ;+ 3 ; 4 ; @todo 5 ; seb 6 ; 7 ;- 8 ;------------------------------------ 1 9 FUNCTION triangule, maskentree, BASIC = basic, COINMONTE = coinmonte, COINDESCEND = coindescend, _extra = ex 2 10 ; -
trunk/SRC/ToBeReviewed/TRIANGULATION/triangule_c.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME:triangule_c 6 ; 7 ; PURPOSE:construit le tableau de triangulation. 8 ; 9 ; L''idee est de 10 ; construire une liste de triangles qui relient les points entre 11 ; eux. Ceci est fait automatiquement avec la fonction TRIANGULATE. 12 ; ICI: 13 ; on tient compte du fait que les points sont disposes sur une grille 14 ; (reguliere ou pas, mais pas destructuree, cad que les points sont 15 ; ecrits suivant une matrice rectangulaire). Un moyen tres simple de 16 ; faire des triangles entre tous les points est alors: 17 ; 18 ; pour chaque point (i,j) de la matrice -sauf ceux de la derniere 19 ; ligne et de la derniere colonne- on on appelle le rectangle 20 ; (i,j) le rectangle forme par les 4 points (i,j), (i+1,j), 21 ; (i,j+1), (i+1,j+1). Pour tracer tous les triangles, il suffit de 22 ; tracer les 2 triangles contenus ds les rectangles (i,j) 23 ; 24 ; au passage on remarque que chaque rectangle (i,j) possede 2 diagonales (si 25 ; si faites un dessin c''est vrai), il y a donc 2 choix possibles pour 26 ; chaque rectangles qd on veut le couper en 2 triangles... 5 ; 6 ; @file_comments 7 ; Construct the triangulation array. 8 ; 9 ; The idea is: construct a list of triangle which link points between them. 10 ; This is automatically done by the function TRIANGULATE 11 ; Here: 12 ; we consider the fact that points are disposed on a grid (regular or not, 13 ; but not unstructured, that is to say that points are written following a 14 ; rectangular matrix). A easy way to do triangles between all points is then: 15 ; 16 ; for each point (i,j) of the matrix -exept those of the last line and of 17 ; the last column- we call rectangle (i,j) the rectangle made of the four 18 ; points (i,j), (i+1,j), (i,j+1), (i+1,j+1). To trace all triangle, we just 19 ; have to trace the 2 triangles contained in rectangles (i,j) 20 ; 21 ; We notice that each rectangle (i,j) have 2 diagonals (it is true... Make a 22 ; drawing to make sure!!), so there are two possible choice for each rectangle 23 ; we want to cut in 2 triangles... 27 24 ; 28 ; C''est grace a ce choix que l''on va pouvoir tracer les cotes avec 29 ; des angles droits. A chaque angle de cote remarquable par 30 ; l''existance d''un unique point terre ou d''un unique point mer sur 31 ; les 4 cotes d''un rectangle (i,j), il faut couper le rectangle 32 ; suivant la diagonale qui qui passe par le point singulier. 25 ; It is thanks to this choice that we will be able to trace coast with right 26 ; angles. At each angle of coast remarkable by the existence of an unique land 27 ; point or of an unique ocean point on one of the four summit of a rectangle (i,j), 28 ; we have to cut the rectangle following the diagonal passing by this point. 33 29 ; 34 ; CATEGORY:pour faire de beaux graphiques masques 35 ; 36 ; CALLING SEQUENCE:res=triangule([mask]) 37 ; 38 ; INPUTS:optionnel:mask c''est le tableau 2d qui sevira a masquer le 39 ; champ que l''on tracera apres avec CONTOUR, 30 ; @categories 31 ; graphic 32 ; 33 ; @param MASKENTREE {in}{optional} 34 ; It is a 2d array which will serve to mask the field we will trace after with CONTOUR, 40 35 ; ...TRIANGULATION=triangule(mask) 41 ; si cet argument n''est pas specifie, la function utilise tmask. 42 ; 43 ; KEYWORD PARAMETERS: 44 ; 45 ; /BASIC: specifie que le masque est sur une grille basice 46 ; (utiliser pour la triangulation ds les coupes verticales et 47 ; des hovmoellers) 48 ; 49 ; /KEEP_CONT: to keep the triangulation even on the continents 50 ; 51 ; COINMONTE=tableau, pour obtenir le tableau de "coins de terre 52 ; montant" a traiter avec completecointerre.pro ds la variable 53 ; tableau plutot que de la faire passer par la variable globale 54 ; twin_corners_up. 55 ; 56 ; COINDESCEND=tableau cf COINMONTE 57 ; 58 ; OUTPUTS: 59 ; res: tableau 2d (3,nbre de triangles). 60 ; chaque ligne de res represente les indices des points 61 ; constituants les sommets d''un triangle. 62 ; cf. comment on trace les triangles ds dessinetri.pro 63 ; 64 ; COMMON BLOCKS: 65 ; common.pro different.pro definetri.pro 66 ; 67 ; SIDE EFFECTS: 68 ; 69 ; RESTRICTIONS:les donnees dont un veut ensuite faire le contour 70 ; doivent etre disposees dans une matrice. Par contre dans la matrice, 71 ; la disposition des points peut ne pas etre irreguliere. Si les 72 ; donnees sont disposees completement de facon irreguliere, utiliser 73 ; TRIANGULE. 74 ; 75 ; EXAMPLE: 76 ; 77 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 36 ; If this argument is not specified, the function use tmask 37 ; 38 ; @keyword BASIC 39 ; Specify that the mask is on a basic grid (use the triangulation for vertical cuts and hovmoellers) 40 ; 41 ; @keyword KEEP_CONT 42 ; To keep the triangulation even on the continents 43 ; 44 ; @keyword COINMONTE 45 ; It is an array. To obtain the array of "ascending land corner" to be treated with 46 ; completecointerre.pro in the variable array instead of make it pass by the global 47 ; variable twin_corners_up. 48 ; 49 ; @keyword COINDESCEND 50 ; It is an array. See COINMONTE 51 ; 52 ; @returns 53 ; res: tableau 2d (3,nbre de triangles). 54 ; Each line of res represent indexes of points constituing summits of a triangle. 55 ; See how we trace triangles in definetri.pro 56 ; 57 ; @uses 58 ; common.pro 59 ; different.pro 60 ; definetri.pro 61 ; 62 ; @restrictions 63 ; Datas whose we want to do the contour must be disposed in a matrix. 64 ; On the other hand, in the matrix, the points's arrangement can not be 65 ; irregular. If it is, use TRIANGULE. 66 ; 67 ; @history 68 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 78 69 ; 26/4/1999 70 ; 71 ; @version 72 ; $Id$ 73 ; 74 ; @todo 75 ; seb L.267->268 je ne pense pas que ce soit ce que tu voulais dire mais 76 ; c'est la traduction de ce qu'il y avait écrit. Correction si besoin. 79 77 ;- 80 78 ;------------------------------------------------------------ … … 83 81 FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, BASIC = basic, KEEP_CONT = keep_cont 84 82 ; 85 86 ; 87 tempsun = systime(1) ; pour key_performance83 compile_opt idl2, strictarrsubs 84 ; 85 tempsun = systime(1) ; For key_performance 88 86 ;--------------------------------------------------------- 89 87 @cm_4mesh 90 88 IF NOT keyword_set(key_forgetold) THEN BEGIN 91 89 @updatenew 92 93 ;------------------------------------------------------------ 94 ; le masque est donne ou il faut prendre tmask?95 ;------------------------------------------------------------ 96 ; 97 98 99 100 101 ; 102 103 ;------------------------------------------------------------ 104 105 106 107 108 109 ;------------------------------------------------------------ 110 ; on va trouver la liste des rectangles (i,j) (reperes par leur coin111 ; en bas a gauche) qu''il faut couper suivant une diagonale descendante112 ; on appellera cette liste: pts_downward90 ENDIF 91 ;------------------------------------------------------------ 92 ; Is the mask given or do we have to take tmask? 93 ;------------------------------------------------------------ 94 ; 95 msk = maskentree 96 taille = size(msk) 97 nx = taille[1] 98 ny = taille[2] 99 ; 100 IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular 101 ;------------------------------------------------------------ 102 if keyword_set(key_periodic)*(nx EQ jpi) $ 103 AND NOT keyword_set(basic) then BEGIN 104 msk = [msk, msk[0, *]] 105 nx = nx+1 106 ENDIF 107 ;------------------------------------------------------------ 108 ; We will find the list of rectangles (i,j)(located by their left 109 ; bottom corner) we have to cut folowing a descendant diagonal. 110 ; We will call this list : pts_downward 113 111 ; 114 115 116 ; on construit le test qui permet de trouver un tel triangle:112 pts_downward = 0 113 114 ; We construct the test which allow to find this triangle : 117 115 ; 118 116 ; … … 124 122 ; msk---------------------shift(msk, -1, 0) 125 123 ; 126 sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;pts qui entourrent le pt en haut a gauche 127 sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;pts qui entourrent le pt en bas a droite 128 129 130 tempdeux = systime(1) ; pour key_performance =2 131 ; pt terre en haut a gauche entoure de pts mer 132 liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 ) 133 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 134 ; pt mer en haut a gauche entoure de pts terre 135 liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1) 136 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 137 ; pt terre en bas a droite entoure de pts mer 138 liste = where( (4-sum2)*(1-shift(msk, -1, 0)) EQ 1) 139 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 140 ; pt mer en bas a droite entoure de pts terre 141 liste = where( (1-sum2)*shift(msk, -1, 0) EQ 1) 142 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 143 undefine, liste 144 ; 145 IF testvar(var = key_performance) EQ 2 THEN $ 146 print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux 147 ; 148 if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 149 tempdeux = systime(1) ; pour key_performance =2 150 ;2 points terre en diagonale montante avec 2 points mer sur la diagonale descendante 151 coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $ 152 *(shift(msk, 0, -1)*shift(msk, -1, 0) EQ 1) ) 153 if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont] 154 ; 155 IF testvar(var = key_performance) EQ 2 THEN $ 156 print, 'temps triangule: trouver coinmont', systime(1)-tempdeux 157 tempdeux = systime(1) ; pour key_performance =2 158 ; 159 ;2 points terre en diagonale descendante avec 2 points mer sur la diagonale montante 160 coindesc = where( ((1-shift(msk, 0, -1))*(1-shift(msk, -1, 0)) $ 161 *msk*shift(msk, -1, -1) EQ 1) ) 162 ; 163 IF testvar(var = key_performance) EQ 2 THEN $ 164 print, 'temps triangule: trouver coindesc', systime(1)-tempdeux 165 ; 166 ENDIF 167 ; 168 if n_elements(pts_downward) EQ 1 then BEGIN 169 tempdeux = systime(1) ; pour key_performance =2 170 ; 171 triang = definetri(nx, ny) 172 ; 173 IF testvar(var = key_performance) EQ 2 THEN $ 174 print, 'temps triangule: definetri', systime(1)-tempdeux 175 coinmont = -1 176 coindesc = -1 177 ENDIF ELSE BEGIN 178 tempdeux = systime(1) ; pour key_performance =2 179 pts_downward = pts_downward[1:n_elements(pts_downward)-1] 180 pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))] 181 ; aucun rectangle ne peut avoir comme coin en bas a gauche un element 182 ; de la derniere colonne ou de la derniere ligne. 183 ; il faut donc enlever ces points si ils ont ete selectionnes dans 184 ; pts_downward. 185 derniere_colonne = (lindgen(ny)+1)*nx-1 186 derniere_ligne = lindgen(nx)+(ny-1)*nx 187 pts_downward =different(pts_downward,derniere_colonne ) 188 pts_downward =different(pts_downward,derniere_ligne ) 189 if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 190 if coinmont[0] NE -1 then begin 124 sum1 = msk+shift(msk, -1, 0)+shift(msk, -1, -1) ;points which surround the left top point. 125 sum2 = msk+shift(msk, 0, -1)+shift(msk, -1, -1) ;points which surround the right bottom point. 126 127 128 tempdeux = systime(1) ; For key_performance =2 129 ; The left top land point surrounded by ocean points 130 liste = where( (4-sum1)*(1-shift(msk, 0, -1)) EQ 1 ) 131 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 132 ; The left top ocean point surrounded by land points 133 liste = where( (1-sum1)*shift(msk, 0, -1) EQ 1) 134 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 135 ; The right bottom land point surrounded by ocean points 136 liste = where( (4-sum2)*(1-shift(msk, -1, 0)) EQ 1) 137 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 138 ; The right bottom ocean point surrounded by land points 139 liste = where( (1-sum2)*shift(msk, -1, 0) EQ 1) 140 if liste[0] NE -1 THEN pts_downward = [pts_downward,liste ] 141 undefine, liste 142 ; 143 IF testvar(var = key_performance) EQ 2 THEN $ 144 print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux 145 ; 146 if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 147 tempdeux = systime(1) ; For key_performance =2 148 ;2 land points in ascendant diagonal with 2 ocean points in descendant diagonal. 149 coinmont = where( (1-msk)*(1-shift(msk, -1, -1)) $ 150 *(shift(msk, 0, -1)*shift(msk, -1, 0) EQ 1) ) 151 if coinmont[0] NE -1 THEN pts_downward = [pts_downward, coinmont] 152 ; 153 IF testvar(var = key_performance) EQ 2 THEN $ 154 print, 'temps triangule: trouver coinmont', systime(1)-tempdeux 155 tempdeux = systime(1) ; pour key_performance =2 156 ; 157 coindesc = where( ((1-shift(msk, 0, -1))*(1-shift(msk, -1, 0)) $ 158 *msk*shift(msk, -1, -1) EQ 1) ) 159 ; 160 ;2 land points in descendant diagonal with 2 ocean points in ascendant diagonal. 161 IF testvar(var = key_performance) EQ 2 THEN $ 162 print, 'temps triangule: trouver coindesc', systime(1)-tempdeux 163 ; 164 ENDIF 165 ; 166 if n_elements(pts_downward) EQ 1 then BEGIN 167 tempdeux = systime(1) ; For key_performance =2 168 ; 169 triang = definetri(nx, ny) 170 ; 171 IF testvar(var = key_performance) EQ 2 THEN $ 172 print, 'temps triangule: definetri', systime(1)-tempdeux 173 coinmont = -1 174 coindesc = -1 175 ENDIF ELSE BEGIN 176 tempdeux = systime(1) ; For key_performance =2 177 pts_downward = pts_downward[1:n_elements(pts_downward)-1] 178 pts_downward = pts_downward[uniq(pts_downward, sort(pts_downward))] 179 ; None rectangle can have an element of the last column or of the 180 ; last line as left bottom corner. 181 ; so we have to remove these points if they has been selected in pts_downward. 182 derniere_colonne = (lindgen(ny)+1)*nx-1 183 derniere_ligne = lindgen(nx)+(ny-1)*nx 184 pts_downward =different(pts_downward,derniere_colonne ) 185 pts_downward =different(pts_downward,derniere_ligne ) 186 if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 187 if coinmont[0] NE -1 then begin 191 188 coinmont =different(coinmont,derniere_colonne ) 192 189 coinmont =different(coinmont,derniere_ligne ) 193 194 190 endif 191 if coindesc[0] NE -1 then begin 195 192 coindesc =different(coindesc,derniere_colonne ) 196 193 coindesc =different(coindesc,derniere_ligne ) 197 endif 198 ENDIF ELSE BEGIN 199 coinmont = -1 200 coindesc = -1 201 ENDELSE 202 IF testvar(var = key_performance) EQ 2 THEN $ 203 print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux 204 ; 205 tempdeux = systime(1) ; pour key_performance =2 206 if pts_downward[0] EQ -1 then triang = definetri(nx, ny) $ 207 ELSE triang = definetri(nx, ny, pts_downward) 208 IF testvar(var = key_performance) EQ 2 THEN $ 209 print, 'temps triangule: definetri', systime(1)-tempdeux 210 ENDELSE 211 ;------------------------------------------------------------ 212 ; on vire les triangles qui ne contiennent que des points terre 213 ;------------------------------------------------------------ 214 ; 215 ; tres bonne idee qui ne marche pas encore a 200% avec IDL 5.2 216 ; ca devrait aller mieux dans les prochaines versions d''IDL... 217 ; 218 if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin 219 tempdeux = systime(1) ; pour key_performance =2 220 ; on enleve les rectangles qui sont entierement dans la terre 221 recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 222 IF testvar(var = key_performance) EQ 2 THEN $ 223 print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 224 225 ; en attendant une version qui marche parfaitement, on est contraint 226 ; de faire un nouveau tri: 227 ; il ne faut pas enlever les rectangles qui n''ont qu''un sommet en 228 ; commun. 194 endif 195 ENDIF ELSE BEGIN 196 coinmont = -1 197 coindesc = -1 198 ENDELSE 199 IF testvar(var = key_performance) EQ 2 THEN $ 200 print, 'temps triangule: menage ds pts_downward coinmont et coindesc', systime(1)-tempdeux 201 ; 202 tempdeux = systime(1) ; For key_performance =2 203 if pts_downward[0] EQ -1 then triang = definetri(nx, ny) $ 204 ELSE triang = definetri(nx, ny, pts_downward) 205 IF testvar(var = key_performance) EQ 2 THEN $ 206 print, 'temps triangule: definetri', systime(1)-tempdeux 207 ENDELSE 208 ;------------------------------------------------------------ 209 ; We delete land points which only contain land points. 210 ;------------------------------------------------------------ 211 ; 212 ; 213 if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin 214 tempdeux = systime(1) ; For key_performance =2 215 ; We delete rectangles which are entirely in the land. 216 recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 217 IF testvar(var = key_performance) EQ 2 THEN $ 218 print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 219 220 ; We do an other sort : 221 ; We have to do not remove rectangles which only have one common summit. 229 222 ; t1 = systime(1) 230 231 232 233 234 tempdeux = systime(1) ; pour key_performance =2235 236 ;if NOT keyword_set(basic) then begin237 vire1 = 0238 vire2 = 0239 while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin240 ; vire sont les rectangles qu''il faut retirer de recsterre (en fait241 ; qu''il faut garder bien qu''ils soient entirement dans la terre)242 vire1 = where( (indice*shift(indice, -1, -1) $243 *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1)244 if vire1[0] NE -1 THEN BEGIN245 indice[vire1] = 0223 indice = intarr(nx, ny) 224 trimask = intarr(nx, ny) 225 trimask[0:nx-2, 0:ny-2] = 1 226 IF recdsterre[0] NE -1 then BEGIN 227 tempdeux = systime(1) ; For key_performance =2 228 indice[recdsterre] = 1 229 if NOT keyword_set(basic) then begin 230 vire1 = 0 231 vire2 = 0 232 while (vire1[0] NE -1 OR vire2[0] NE -1) ne 0 do begin 233 ; Delete rectangles we have to remove from recsterre (in fact those we have 234 ; to keep although they ar eentirely in the land). 235 vire1 = where( (indice*shift(indice, -1, -1) $ 236 *(1-shift(indice, 0, -1))*(1-shift(indice, -1, 0))*trimask) EQ 1) 237 if vire1[0] NE -1 THEN BEGIN 238 indice[vire1] = 0 246 239 ; indice[vire1+nx+1] = 0 247 endif248 249 vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $250 *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1)251 if vire2[0] NE -1 THEN BEGIN252 indice[vire2+1] = 0240 endif 241 242 vire2 = where( ((1-indice)*(1-shift(indice, -1, -1)) $ 243 *shift(indice, 0, -1)*shift(indice, -1, 0)*trimask) EQ 1) 244 if vire2[0] NE -1 THEN BEGIN 245 indice[vire2+1] = 0 253 246 ; indice[vire2+nx] = 0 254 endif255 endwhile256 IF testvar(var = key_performance) EQ 2 THEN $257 print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux258 ;endif259 indice[*, ny-1] = 1 ; la deriere colonne te la derniere ligne260 indice[nx-1, *] = 1 ; ne peuvent definir derectangle.261 ; 262 tempdeux = systime(1) ; pour key_performance =2263 264 ; on recupere les numeros des triangles que l'' on va garder265 266 267 247 endif 248 endwhile 249 IF testvar(var = key_performance) EQ 2 THEN $ 250 print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux 251 endif 252 indice[*, ny-1] = 1 ; The last column and the last line 253 indice[nx-1, *] = 1 ; can not define any rectangle. 254 ; 255 tempdeux = systime(1) ; For key_performance =2 256 recgarde = where(indice EQ 0) 257 ; We recuperate numbers of triangles we will keep. 258 trigarde = 2*[recgarde-recgarde/nx] 259 trigarde = transpose(temporary(trigarde)) 260 trigarde = [trigarde, trigarde+1] 268 261 ; 269 270 262 triang = triang[*, temporary(trigarde[*])] 263 IF testvar(var = key_performance) EQ 2 THEN $ 271 264 print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux 272 273 265 endif 266 endif 274 267 ; print, 'temps tri triangles', systime(1)-t1 275 268 ;------------------------------------------------------------ 276 ; quand key_periodic eq 1, triang est une liste d''indice d'un277 ; tableau qui a une colonne de trop.278 ; il faut ramener ca a la matrice initiale en mettant les indivces de279 ; la derniere colonne egaux a ceux de la derniere colonne...280 ;------------------------------------------------------------ 281 tempdeux = systime(1) ; pour key_performance =2282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 ;------------------------------------------------------------ 314 315 316 ;------------------------------------------------------------ 317 269 ; When key_periodic equal 1, triang is a list of indexes's array which 270 ; have a surplus column. 271 ; We have to put it back to the initial matrix by putting indexes of 272 ; the last column equal to these of the last column... 273 ;------------------------------------------------------------ 274 tempdeux = systime(1) ; For key_performance =2 275 if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 276 AND NOT keyword_set(basic) then BEGIN 277 indicey = triang/nx 278 indicex = triang-indicey*nx 279 nx = nx-1 280 liste = where(indicex EQ nx) 281 if liste[0] NE -1 then indicex[liste] = 0 282 triang = indicex+nx*indicey 283 nx = nx+1 284 if coinmont[0] NE -1 then begin 285 indicey = coinmont/nx 286 indicex = coinmont-indicey*nx 287 nx = nx-1 288 liste = where(indicex EQ nx) 289 if liste[0] NE -1 THEN indicex[liste] = 0 290 coinmont = indicex+nx*indicey 291 nx = nx+1 292 endif 293 if coindesc[0] NE -1 then begin 294 indicey = coindesc/nx 295 indicex = coindesc-indicey*nx 296 nx = nx-1 297 liste = where(indicex EQ nx) 298 if liste[0] NE -1 THEN indicex[liste] = 0 299 coindesc = indicex+nx*indicey 300 nx = nx+1 301 endif 302 endif 303 IF testvar(var = key_performance) EQ 2 THEN $ 304 print, 'temps triangule: finitions', systime(1)-tempdeux 305 306 ;------------------------------------------------------------ 307 if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 308 if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 309 ;------------------------------------------------------------ 310 IF NOT keyword_set(key_forgetold) THEN BEGIN 318 311 @updateold 319 320 321 322 323 312 ENDIF 313 314 IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun 315 316 return, triang 324 317 325 318 END -
trunk/SRC/ToBeReviewed/TRIANGULATION/triangule_e.pro
r134 r150 3 3 ;------------------------------------------------------------ 4 4 ;+ 5 ; NAME:triangule_e 6 ; 7 ; PURPOSE:buid the triangulation for a E-grid type 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) 5 ; @file_comments 6 ; Buid the triangulation for a E-grid type 7 ; 8 ; @categories 9 ; graphic 10 ; 11 ; @param MASKENTREE {in}{optional} 12 ; It is a 2d array which will serve to mask the field we will trace after with CONTOUR, 13 ; ...TRIANGULATION=triangule(mask) 14 ; If this argument is not specified, the function use tmask 15 ; 16 ; @keyword BASIC 17 ; Specify that the mask is on a basic grid (use the triangulation for vertical cuts and hovmoellers) 18 ; 19 ; @keyword COINMONTE 20 ; It is an array. To obtain the array of "ascending land corner" to be treated with 21 ; completecointerre.pro in the variable array instead of make it pass by the global 22 ; variable twin_corners_up. 23 ; 24 ; @keyword COINDESCEND 25 ; It is an array. See COINMONTE 26 ; 27 ; @keyword SHIFTED 28 ; 29 ; @uses 30 ; common.pro 31 ; 32 ; @history 33 ; Sebastien Masson (smasson@lodyc.jussieu.fr) 28 34 ; june 2001 35 ; 36 ; @version 37 ; $Id$ 38 ; 39 ; @todo 40 ; seb L.152->153 je ne pense pas que ce soit ce que tu voulais dire mais 41 ; c'est la traduction de ce qu'il y avait écrit. Correction si besoin. 29 42 ;- 30 43 ;------------------------------------------------------------ … … 42 55 ENDIF 43 56 ;--------------------------------------------------------- 44 tempsun = systime(1) ; pour key_performance45 ;------------------------------------------------------------ 46 ; le masque est donne ou il faut prendre tmask?57 tempsun = systime(1) ; For key_performance 58 ;------------------------------------------------------------ 59 ; Is the mask given or do we have to take tmask? 47 60 ;------------------------------------------------------------ 48 61 ; … … 138 151 ; ; 139 152 ;------------------------------------------------------------ 140 ; quand key_periodic eq 1, triang est une liste d''indice d'un141 ; tableau qui a une colonne de trop.142 ; il faut ramener ca a la matrice initiale en mettant les indivces de143 ; la derniere colonne egaux a ceux de la derniere colonne...144 ;------------------------------------------------------------ 145 tempdeux = systime(1) ; pour key_performance =2153 ; When key_periodic equal 1, triang is a list of indexes's array which 154 ; have a surplus column. 155 ; We have to put it back to the initial matrix by putting indexes of 156 ; the last column equal to these of the last column... 157 ;------------------------------------------------------------ 158 tempdeux = systime(1) ; For key_performance =2 146 159 if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 147 160 AND NOT keyword_set(basic) then BEGIN … … 176 189 177 190 ;------------------------------------------------------------ 178 ;if arg_present(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont179 ;if arg_present(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc180 ; 181 ;IF NOT keyword_set(key_forgetold) THEN BEGIN182 ;@updateold183 ;ENDIF191 if arg_present(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 192 if arg_present(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 193 194 IF NOT keyword_set(key_forgetold) THEN BEGIN 195 @updateold 196 ENDIF 184 197 ;------------------------------------------------------------ 185 198
Note: See TracChangeset
for help on using the changeset viewer.