; $Id: map_proj_image.pro,v 1.4 2004/07/08 23:45:53 chris Exp $ ; ; Copyright (c) 2004-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;------------------------------------------------------------------------- ;+ ; NAME: ; map_proj_image ; ; PURPOSE: ; This function returns the Image_Orig image warped to fit ; the current map. Image_Orig must be centered at 0. This ; routine works in image space. ; ; SYNTAX: ; Result = MAP_PROJ_IMAGE(Image_Orig [, Range] ) ; ; INPUT: ; Image_Orig: A two-dimensional array representing geographical ; image data. ; ; Range: An optional 4-element vector giving the image extents as: ; [Lonmin, Latmin, Lonmax, Latmax]. ; ; KEYWORDS: ; BILINEAR: Set this keyword to use bilinear interpolation. The ; default is nearest neighbor. ; ; DIMENSIONS: Set this keyword to a two-element vector giving ; the width and height of the desired Result. ; The UVRANGE keyword may be used to retrieve ; the resulting range of the UV coordinates. ; If DIMENSIONS is not specified then the Result has the ; same dimensions as the input image. ; ; IMAGE_STRUCTURE: Set this keyword to a !MAP structure as returned ; from MAP_PROJ_INIT, giving the map projection of the ; original image. If this keyword is present, the original image ; is assumed to be in meters. In this case, the image is first ; warped from the IMAGE_STRUCTURE projection to degrees lat/lon, ; and then wraped to the final MAP_STRUCTURE. ; If this keyword is omitted, the original image is assumed to be ; in degrees longitude/latitude (geographic coordinates). ; ; MAP_STRUCTURE: Set this keyword to a !MAP structure as returned ; from MAP_PROJ_INIT, giving the map projection to which to ; warp the image. If this keyword is omitted, ; the !MAP system variable is used. ; ; MASK: Set this keyword to a named variable in which to return ; a byte array of the same dimensions as the Result, containing ; a mask of the good values. Values in the Result that were ; set to MISSING will have a mask value of zero, while all ; other mask values will be one. ; ; MISSING: Set this keyword to the value to set areas outside ; the valid map coordinates. If omitted, areas outside the map ; are set to 0. ; ; MAX_VALUE = values in Image_Orig greater than or equal to MAX_VALUE ; are considered missing. Pixels in the output image ; that depend upon missing pixels will be set to MISSING. ; ; MIN_VALUE = values in Image_Orig less than or equal to MIN_VALUE ; are considered missing. ; ; UVRANGE: Set this keyword to a named variable in which to return ; the range of the UV coordinates as: [umin, vmin, umax, vmax]. ; ; OUTPUT: ; The warped image is returned. ; ; PROCEDURE: An image space algorithm is used, so the time required ; is roughly proportional to the size of the final image. ; For each pixel in the Result, the inverse coordinate transform is ; applied to obtain lat/lon. The lat/lon coordinates are then scaled ; into image pixel coordinates, and these coordinates are then ; interpolated from Image values. ; ; MODIFICATION HISTORY: ; Written: CT, RSI, April 2004. ; ;- ;------------------------------------------------------------------------- ; Given a map structure and a latlon range, determine the corresponding ; range of UV coordinates. ; function map_proj_image_uvsize, xminIn, xmaxIn, yminIn, ymaxIn, $ IMAGE_STRUCTURE=imageStructure, $ MAP_STRUCTURE=mapStruct compile_opt idl2, hidden hasImageStruct = N_TAGS(imageStructure) gt 0 ; Subdivisions across lat/lon range. nx = 61 ny = 31 if (hasImageStruct) then begin xmin = xminIn xmax = xmaxIn ymin = yminIn ymax = ymaxIn dx = (xmax-xmin)/(nx-1) dy = (ymax-ymin)/(ny-1) xcoord = REFORM(REBIN(DINDGEN(nx)*dx + xmin, nx, ny), nx*ny) ycoord = REFORM(REBIN(DINDGEN(1,ny)*dy + ymin, nx, ny), nx*ny) ; First convert from original UV to lonlat. lon = MAP_PROJ_INVERSE(xcoord, ycoord, MAP_STRUCTURE=imageStructure) lat = REFORM(lon[1,*]) lon = REFORM(lon[0,*]) ; Filter out NaNs. good = WHERE(FINITE(lon) and FINITE(lat)) lon = lon[good] lat = lat[good] endif else begin xmin = xminIn > mapStruct.ll_box[1] xmax = xmaxIn < mapStruct.ll_box[3] ymin = yminIn > mapStruct.ll_box[0] ymax = ymaxIn < mapStruct.ll_box[2] dx = (xmax-xmin)/(nx-1) dy = (ymax-ymin)/(ny-1) lon = REFORM(REBIN(DINDGEN(nx)*dx + xmin, nx, ny), nx*ny) lat = REFORM(REBIN(DINDGEN(1,ny)*dy + ymin, nx, ny), nx*ny) endelse ; Convert from lonlat to UV. u = MAP_PROJ_FORWARD(lon, lat, MAP_STRUCTURE=mapStruct) v = REFORM(u[1,*]) u = REFORM(u[0,*]) ; Filter out NaNs. good = WHERE(FINITE(u) and FINITE(v), ngood) if (~ngood) then $ return, [0d, 0, 1, 1] u = u[good] v = v[good] minu = MIN(u, MAX=maxu) minv = MIN(v, MAX=maxv) return, [minu, minv, maxu, maxv] end ;------------------------------------------------------------------------- ; Return an array of 1's where the data are outside the range of min_value ; to max_value. Max_value and/or min_value may be undefined. If both are ; undefined, return a -1. ; function map_proj_image_missing, image_orig, max_value, min_value compile_opt idl2, hidden if (N_ELEMENTS(max_value) eq 1) then begin if (N_ELEMENTS(min_value) eq 1) then $ return, (Image_orig ge max_value) or (Image_orig le min_value) return, Image_orig ge max_value endif if (N_ELEMENTS(min_value) eq 1) then $ return, Image_orig le min_value return, -1 end ;------------------------------------------------------------------------- function map_proj_image, Image_Orig, rangeIn, $ BILINEAR = bilin, $ DIMENSIONS=dimensions, $ IMAGE_STRUCTURE=imageStructure, $ MAP_STRUCTURE=mapStruct, $ MASK=mask, $ MAX_VALUE = max_value, $ MIN_VALUE=min_value, $ MISSING = missingIn, $ UVRANGE=uvrange, $ XINDEX=xindex, $ ; currently undocumented YINDEX=yindex ; currently undocumented compile_opt idl2 ON_ERROR,2 dims = SIZE(Image_Orig, /DIMENSIONS) if (N_ELEMENTS(dims) ne 2) then $ MESSAGE, "Image must be a two- dimensional array." s1 = dims[0] ; # of columns s2 = dims[1] ; # of rows if (s1 le 1) || (s2 le 1) THEN $ message, 'Each dimension must be greater than 1." ; Missing data value should equal the background or user-supplied value. missing = (N_ELEMENTS(missingIn) eq 1) ? missingIn : 0 if (N_ELEMENTS(dimensions) eq 2) then begin nx = LONG(dimensions[0]) ny = LONG(dimensions[1]) endif else begin nx = s1 ny = s2 endelse if (N_TAGS(mapStruct) eq 0) then $ mapStruct = !MAP hasImageStruct = N_TAGS(imageStructure) gt 0 if (N_ELEMENTS(rangeIn) eq 4) then begin xmin = DOUBLE(rangeIn[0]) ymin = DOUBLE(rangeIn[1]) xmax = DOUBLE(rangeIn[2]) ymax = DOUBLE(rangeIn[3]) endif else begin if (hasImageStruct) then $ MESSAGE, 'Range must be specified with IMAGE_STRUCTURE keyword.' ; Assume image covers entire globe. xmin = -180d ymin = -90d xmax = 180d ymax = 90d endelse ; Scale from lat/lon to pixels sx = ((s1-1.)/(xmax - xmin)) sy = ((s2-1.)/(ymax - ymin)) ; Does image wrap around globe? if (~hasImageStruct) then begin wrap = ((xmin - 360./s1 + 720.) mod 360.) - ((xmax + 720.) mod 360.) wrap = abs(wrap) lt 1e-4 ;Allow for roundoff endif else wrap = 0 ; If user passed in cached values for XINDEX and YINDEX then use them. if (N_ELEMENTS(xindex) gt 0 && N_ELEMENTS(yindex) gt 0) then $ goto, skipover ; Find the extent of the our limits in the map on the screen by ; making a n x n array of lon/lats spaced over the extent of ; the image and saving the extrema. uvrange = MAP_PROJ_IMAGE_UVSIZE(xmin, xmax, ymin, ymax, $ IMAGE_STRUCTURE=imageStructure, $ MAP_STRUCTURE=mapStruct) ; U and V Cartesian coordinates. u = DINDGEN(nx)*((uvrange[2] - uvrange[0])/(nx-1)) + uvrange[0] v = DINDGEN(1, ny)*((uvrange[3] - uvrange[1])/(ny-1)) + uvrange[1] uvIn = [REFORM(REBIN(u, nx, ny), 1, nx*ny), $ REFORM(REBIN(v, nx, ny), 1, nx*ny)] ; Convert from UV to lonlat. lonlat = MAP_PROJ_INVERSE(uvIn, MAP_STRUCTURE=mapStruct) ; The GCTP projections tend to collapse invalid UV points ; down to "valid" lonlat. To filter out these bad points, ; do a forward transform of our lonlat, and compare the ; resulting UV values with the originals. ; First filter out NaN's to avoid overflow errors. bad = WHERE(~FINITE(lonlat), nbad) if (nbad gt 0) then $ lonlat[bad] = 0 ; arbitrary value lon = lonlat[0,*] lat = lonlat[1,*] ; Throw away out-of-bounds latitude values. ibad = ABS(lat) gt 90 bad = WHERE(ibad, nbad) if (nbad gt 0) then $ lonlat[1,bad] = 0 ; If the bounding box is less than the image range, ; then clip the lonlats. ibad or= (lon lt mapStruct.ll_box[1]) ibad or= (lon gt mapStruct.ll_box[3]) ibad or= (lat lt mapStruct.ll_box[0]) ibad or= (lat gt mapStruct.ll_box[2]) uvOut = MAP_PROJ_FORWARD(lonlat, MAP_STRUCTURE=mapStruct) ; Calculate relative error between uvIn/Out. diff = TEMPORARY(uvOut) - uvIn ; Avoid division by zero. zero = WHERE(~uvIn, nzero) if (nzero gt 0) then $ uvIn[zero] = 1 diff = ABS(TEMPORARY(diff)/TEMPORARY(uvIn)) ; Somewhat arbitrary cutoff of "bad" values. ibad or= MAX(TEMPORARY(diff) gt 0.1d, DIM=1) ; Mark all bad values. bad = WHERE(ibad, nbad) if (hasImageStruct) then begin xcoord = MAP_PROJ_FORWARD(lonlat, MAP_STRUCTURE=imageStructure) if (nbad gt 0) then $ xcoord[*, bad] = 1d10 ycoord = REFORM(xcoord[1,*], nx, ny) xcoord = REFORM(xcoord[0,*], nx, ny) endif else begin if (nbad gt 0) then $ lonlat[*, bad] = 1d10 lat = reform(lonlat[1,*], nx, ny) ; Separate lat/lon lon = reform(lonlat[0,*], nx, ny) lonlat = 0 ; free memory w = where(lon lt xmin, count) ;Handle longitude wrap-around while count gt 0 do begin lon[w] += 360.0 w = where(lon lt xmin, count) endwhile xcoord = TEMPORARY(lon) ycoord = TEMPORARY(lat) endelse skipover: ; Now interpolate the screen image from the original. if KEYWORD_SET(bilin) then begin ; Didn't get cached values? if (N_ELEMENTS(xindex) eq 0) then begin xindex = (TEMPORARY(xcoord) - xmin) * sx ;To pixels yindex = (TEMPORARY(ycoord) - ymin) * sy endif ; Verify that we indeed have pixels that wrap around. if wrap then begin ; We will use col1 below if necessary. col1 = where(xindex gt (s1-1), nwrap) if (~nwrap) then $ wrap = 0 endif ; Points off the map will be set to missing value. Image_Warp = INTERPOLATE(Image_Orig, xindex, yindex, MISSING = missing) ; If the image wraps around the globe, we must treat those pixels ; after the last column and before the first column specially. if (wrap) then begin threecol = [Image_Orig[s1-1,*], Image_Orig[0:1,*]] col1x = xindex[col1] - (s1-1) ;Interpolate value ; Add in points that wrapped in xindex between s1-1 and s1. Image_Warp[col1] = INTERPOLATE(threecol, col1x, yindex[col1], $ MISSING = missing) endif badb = MAP_PROJ_IMAGE_MISSING(image_orig, max_value, min_value) if badb[0] ne -1 then begin ;Missing data value? ; Interpolate location of missing data values onto new grid, ; then force to missing. mask = ~INTERPOLATE(float(badb), xindex, yindex, MISSING=1) if (wrap) then begin mask[col1] = ~INTERPOLATE(float([badb[s1-1,*], badb[0:1,*]]), $ col1x, yindex[col1], MISSING=1) endif bad = where(~mask, nbad) if (nbad gt 0) then $ Image_Warp[bad] = missing endif else begin ; no badb ; Don't bother with MASK unless user requested it. if ARG_PRESENT(mask) then begin ; Do a fake interpolation of an array of ones to find the points ; off the map. Set these to missing=zero. mask = INTERPOLATE(REPLICATE(1b, s1, s2), xindex, yindex, $ MISS = 0) if (wrap) then begin mask[col1] = INTERPOLATE(REPLICATE(1b, 3, s2), $ col1x, yindex[col1], MISS = 0) endif endif endelse ; no badb endif else begin ; Nearest neighbor ; Didn't get cached values? if (N_ELEMENTS(xindex) eq 0) then begin ; This is the same as: xindex = (lon-xmin) * sx + 0.5, but faster for arrays. xindex = (TEMPORARY(xcoord) - (xmin - 0.5/sx)) * sx yindex = (TEMPORARY(ycoord) - (ymin - 0.5/sy)) * sy if wrap then begin needWrap = where(xindex ge s1 and xindex lt (2*s1), count) if count gt 0 then $ xindex[TEMPORARY(needWrap)] -= s1 endif endif ; It's okay to index using out-of-bounds values (they will be ; clipped), but we need to save the bad locations so we can ; set to missing at the end. mask = (xindex ge 0 and xindex lt s1 and yindex ge 0 and yindex lt s2) Image_Warp = Image_Orig[xindex, yindex] badb = MAP_PROJ_IMAGE_MISSING(Image_Warp, max_value, min_value) if ((badb[0] ne -1) && ~ARRAY_EQUAL(badb, 1)) then begin mask and= ~badb endif bad = where(~mask, offmap) if (offmap gt 0) then $ Image_Warp[bad] = missing endelse ; Nearest neighbor return, Image_Warp end