;+ ; ; @file_comments ; compute the weight and address needed to interpolate data from ; an "irregular 2D grid" (defined as a grid made of quadrilateral cells) ; to any grid using the bilinear method ; ; @categories ; Interpolation ; ; @param olonin {in}{required}{type=2d array} ; longitude of the input data ; ; @param olat {in}{required}{type=2d array} ; latitude of the input data ; ; @param omsk {in}{required}{type=2d array or -1} ; land/sea mask of the input data ; put -1 if input data are not masked ; ; @param alonin {in}{required}{type=2d array} ; longitude of the output data ; ; @param alat {in}{required}{type=2d array} ; latitude of the output data ; ; @param amsk {in}{required}{type=2d array or -1} ; land/sea mask of the output data ; put -1 if output data are not masked ; ; @param weig {out}{type=2d array} ; (see ADDR) ; ; @param addr {out}{type=2d array} ; 2D arrays, weig and addr are the weight and addresses used to ; perform the interpolation: ; dataout = total(weig*datain[addr], 1) ; dataout = reform(dataout, jpia, jpja, /over) ; ; @restrictions ; - the input grid must be an "irregular 2D grid", defined as a grid made ; of quadrilateral cells which corners positions are defined with olonin and olat ; - We supposed the data are located on a sphere, with a periodicity along ; the longitude ; - to perform the bilinear interpolation within quadrilateral cells, we ; first morph the cell into a square cell and then compute the bilinear ; interpolation. ; - if some corners of the cell are land points, their weight is set to 0 ; and the weight is redistributed on the remaining "water" corners ; - points located out of the southern and northern boundaries or in cells ; containing only land points are set the same value as their closest neighbor ; ; @history ; June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr) ; ; @version ; $Id$ ; ;- ; PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk, weig, addr ; compile_opt idl2, strictarrsubs ; jpia = (size(alonin, /dimensions))[0] jpja = (size(alonin, /dimensions))[1] ; jpio = (size(olonin, /dimensions))[0] jpjo = (size(olonin, /dimensions))[1] ; mask check IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN BEGIN omsk = replicate(1b, jpio, jpjo) IF jpio EQ 182 AND jpjo EQ 149 THEN BEGIN lontmp = (olonin[1:180, *] + 3600) MOD 360 lattmp = olat[1:180, *] bad = abs(abs(lontmp - shift(lontmp, 1, 0)) - 180) LT 176 $ OR abs(abs(lontmp - shift(lontmp, 0, 1)) - 180) LT 176 $ OR abs(abs(lontmp - shift(lontmp, 1, 1)) - 180) LT 176 $ OR abs(abs(lattmp - shift(lattmp, 1, 0)) - 180) LT 176 $ OR abs(abs(lattmp - shift(lattmp, 0, 1)) - 180) LT 176 $ OR abs(abs(lattmp - shift(lattmp, 1, 1)) - 180) LT 176 bad = bad AND lattmp LT 50 omsk[1:180, 1:148] = 1b - bad[*, 1:148] omsk[0, *] = omsk[180, *] omsk[181, *] = omsk[1, *] ENDIF ELSE omsk = replicate(1b, jpio, jpjo) ENDIF IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) IF n_elements(omsk) NE jpio*jpjo THEN BEGIN ras = report('input grid mask do not have the good size') stop ENDIF IF n_elements(amsk) NE jpia*jpja THEN BEGIN ras = report('output grid mask do not have the good size') stop ENDIF ; ; longitude, between 0 and 360 alon = alonin MOD 360 out = where(alon LT 0) IF out[0] NE -1 THEN alon[out] = alon[out]+360 olon = olonin MOD 360 out = where(olon LT 0) IF out[0] NE -1 THEN olon[out] = olon[out]+360 ; ; we work only on the ouput grid water points awater = where(amsk EQ 1) nawater = n_elements(awater) ; ; define all cells that have corners located at olon, olat ; we define the cell with the address of all corners ; ; 3 2 ; +------+ ; | | ; | | ; | | ; +------+ ; 0 1 ; alladdr = lindgen(jpio, jpjo-1) celladdr = lonarr(4, jpio*(jpjo-1)) celladdr[0, *] = alladdr celladdr[1, *] = shift(alladdr, -1) celladdr[2, *] = shift(alladdr + jpio, -1) celladdr[3, *] = alladdr + jpio ; ; list the cells that have at least 1 ocean point as corner good = where(total(omsk[celladdr], 1) GT 0) ; keep only those cells celladdr = celladdr[*, temporary(good)] ; xcell = olon[celladdr] minxcell = min(xcell, dimension = 1, max = maxxcell) ycell = olat[celladdr] minycell = min(ycell, dimension = 1, max = maxycell) ; foraddr: address of the ocean water cell associated to each atmosphere water point foraddr = lonarr(nawater) ; forweight: x/y position of the atmosphere water point in the ocean water cell forweight = dblarr(nawater, 2) ; ; Loop on all the water point of the atmosphere ; We look for which ocean water cell contains the atmosphere water point ; delta = max([(360./jpio), (180./jpjo)])* 4. FOR n = 0L, nawater-1 DO BEGIN ; control print IF (n MOD 5000) EQ 0 THEN print, n ; longitude and latitude of the atmosphere water point xx = alon[awater[n]] yy = alat[awater[n]] ; 1) we reduce the number of ocean cells for which we will try to know if ; xx,yy is inside. CASE 1 OF ; if we are near the norh pole yy GE (90-delta):BEGIN lat1 = 90-2*delta good = where(maxycell GE lat1) onsphere = 1 END ; if we are near the longitude periodicity area xx LE delta OR xx GE (360-delta):BEGIN lat1 = yy-delta lat2 = yy+delta good = where((minxcell LE 2*delta OR maxxcell GE (360-2*delta)) AND maxycell GE lat1 AND minycell LE lat2) onsphere = 1 END ; other cases ELSE:BEGIN lon1 = xx-delta lon2 = xx+delta lat1 = yy-delta lat2 = yy+delta good = where(maxxcell GE lon1 AND minxcell LE lon2 AND maxycell GE lat1 AND minycell le lat2) ; ORCA cases : orca grid is irregular only northward of 40N CASE 1 OF jpio EQ 92 AND (jpjo EQ 76 OR jpjo EQ 75 OR jpjo EQ 74 ):onsphere = yy GT 40 jpio EQ 180 AND (jpjo EQ 149 OR jpjo EQ 148 OR jpjo EQ 147 ):onsphere = yy GT 40 jpio EQ 720 AND (jpjo EQ 522 OR jpjo EQ 521 OR jpjo EQ 520 ):onsphere = yy GT 40 jpio EQ 1440 AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 ELSE:onsphere = 1 ENDCASE END ENDCASE ; we found a short list of possible ocean water cells containing the atmosphere water point IF good[0] NE -1 THEN BEGIN ; in which cell is located the atmosphere water point? ; Warning! inquad use clockwise quadrilateral definition ind = inquad(xx, yy $ , xcell[0, good], ycell[0, good] $ , xcell[3, good], ycell[3, good] $ , xcell[2, good], ycell[2, good] $ , xcell[1, good], ycell[1, good] $ , onsphere = onsphere, newcoord = newcoord, /noprint) ; keep only the first cell (if the atmospheric point was located in several ocean cells) ind = ind[0] ; we found one ocean water cell containing the atmosphere water point IF ind NE -1 THEN BEGIN ind = good[ind] ; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1) ; in addition we get the coordinates of the atmospheric point in this new morphed square IF onsphere THEN BEGIN ; Warning! quadrilateral2square use anticlockwise quadrilateral definition xy = quadrilateral2square(newcoord[0, 0], newcoord[1, 0] $ , newcoord[0, 3], newcoord[1, 3] $ , newcoord[0, 2], newcoord[1, 2] $ , newcoord[0, 1], newcoord[1, 1] $ , newcoord[0, 4], newcoord[1, 4]) ENDIF ELSE BEGIN xy = quadrilateral2square(xcell[0, ind], ycell[0, ind] $ , xcell[1, ind], ycell[1, ind] $ , xcell[2, ind], ycell[2, ind] $ , xcell[3, ind], ycell[3, ind], xx, yy) ENDELSE ; take care of rounding errors... zero = where(abs(xy) LT 1e-4) IF zero[0] NE -1 THEN xy[zero] = 0 one = where(abs(1-xy) LT 1e-4) IF one[0] NE -1 THEN xy[one] = 1 ; some checks... tmpmsk = omsk[celladdr[*, ind]] CASE 1 OF xy[0] LT 0 OR xy[0] GT 1 : stop xy[1] LT 0 OR xy[1] GT 1 : stop xy[0] EQ 0 AND xy[1] EQ 0 AND tmpmsk[0] EQ 0 : foraddr[n] = -1 xy[0] EQ 1 AND xy[1] EQ 0 AND tmpmsk[1] EQ 0 : foraddr[n] = -1 xy[0] EQ 1 AND xy[1] EQ 1 AND tmpmsk[2] EQ 0 : foraddr[n] = -1 xy[0] EQ 0 AND xy[1] EQ 1 AND tmpmsk[3] EQ 0 : foraddr[n] = -1 xy[0] EQ 0 AND (tmpmsk[0]+tmpmsk[3]) EQ 0 : foraddr[n] = -1 xy[0] EQ 1 AND (tmpmsk[1]+tmpmsk[2]) EQ 0 : foraddr[n] = -1 xy[1] EQ 0 AND (tmpmsk[0]+tmpmsk[1]) EQ 0 : foraddr[n] = -1 xy[1] EQ 1 AND (tmpmsk[2]+tmpmsk[3]) EQ 0 : foraddr[n] = -1 ELSE: BEGIN ; we keep its address foraddr[n] = ind ; keep the new coordinates forweight[n, 0] = xy[0] forweight[n, 1] = xy[1] END ENDCASE ENDIF ELSE foraddr[n] = -1 ENDIF ELSE foraddr[n] = -1 ENDFOR ; do we have some water atmospheric points that are not located in an water oceanic cell? bad = where(foraddr EQ -1) IF bad[0] NE -1 THEN BEGIN ; yes? ; we look for neighbor water atmospheric point located in water oceanic cell badaddr = awater[bad] good = where(foraddr NE -1) ; list the atmospheric points located in water oceanic cell goodaddr = awater[good] ; there longitude and latitude goodlon = alon[goodaddr] goodlat = alat[goodaddr] ; for all the bad points, look for a neighbor neig = lonarr(n_elements(bad)) FOR i = 0, n_elements(bad)-1 DO BEGIN neig[i] = (neighbor(alon[badaddr[i]], alat[badaddr[i]], goodlon, goodlat, /sphere))[0] ENDFOR ; get the address regarding foraddr neig = good[neig] ; associate each bad point with its neighbor (get its address and weight) foraddr[bad] = foraddr[neig] forweight[bad, *] = forweight[neig, *] endif ; transform the address of the ocean cell into the address of its 4 corners newaaddr = celladdr[*, temporary(foraddr)] ; now we compute the weight to give at each corner newaweig = dblarr(4, nawater) a = reform(forweight[*, 0], 1, nawater) b = reform(forweight[*, 1], 1, nawater) forweight = -1 ; free memory newaweig = [(1-a)*(1-b), (1-b)*a, a*b, (1-a)*b] a = -1 & b = -1 ; free memory ; mask the weight to suppress the corner located on land newaweig = newaweig*((omsk)[newaaddr]) totalweig = total(newaweig, 1) ; for cell with some land corner, ; we have to redistribute the weight on the reaining water corners ; weights normalization totalweig = total(newaweig, 1) IF min(totalweig, max = ma) EQ 0 then stop IF ma GT 1 then stop newaweig = newaweig/(replicate(1., 4)#totalweig) totalweig = total(newaweig, 1) ; weights weig = dblarr(4, jpia*jpja) weig[*, awater] = temporary(newaweig) ; address addr = dblarr(4, jpia*jpja) addr[*, awater] = temporary(newaaddr) ; RETURN END