; docformat = 'rst' ; ; PURPOSE: ; The purpose of this function is to warp an image into a map projection, given ; latitude and longitude values for each data point. It is similar to MAP_PATCH in ; IDL. ; ;******************************************************************************************; ; ; ; Copyright (c) 2012, by Fanning Software Consulting, Inc. All rights reserved. ; ; ; ; Redistribution and use in source and binary forms, with or without ; ; modification, are permitted provided that the following conditions are met: ; ; ; ; * Redistributions of source code must retain the above copyright ; ; notice, this list of conditions and the following disclaimer. ; ; * Redistributions in binary form must reproduce the above copyright ; ; notice, this list of conditions and the following disclaimer in the ; ; documentation and/or other materials provided with the distribution. ; ; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ; ; contributors may be used to endorse or promote products derived from this ; ; software without specific prior written permission. ; ; ; ; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ; ; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ; ; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ; ; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ; ; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ; ; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ; ; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ; ; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ; ; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ; ; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ; ;******************************************************************************************; ; ;+ ; The purpose of this function is to warp an image into a map projection, given ; latitude and longitude values for each data point. It is the Coyote Graphics ; equivalent of MAP_PATCH in IDL. ; ; Please note: I have only used and tested this function with the well-behaved ; data I am using in my own research and for which I needed this functionality. ; I have no doubt there might be problems with less well-behaved data sets. ; If you discover a problem with your own data, please let me know and I'll take ; another look at this function. DWF. ; ; :Categories: ; Graphics, Map Projections ; ; :Returns: ; An output 2D grid in which the data points have been warped into the ; particular map projection at the requested pixel resolution. ; ; :Params: ; data: in, required, type=numerical ; A vector or 2D data set. Longitude and latitude values must be present ; (or easily calculated) for each element in the data array. ; lons: in, required, type=float ; If data is a vector, a vector of longitude values corresponding to each ; data value. Values must be in the range -180 to 360. If data is 2D, either ; a 1D vector or a 2D array of corresponding longitude values. If data is 2D, ; and the LONS parameter is missing, a vector of appropriate length scaled into ; the range -180.0 to 180.0 will be created. ; lats: in, required, type=float ; If data is a vector, a vector of latitude values corresponding to each ; data value. Values must be in the range -90 to 90. If data is 2D, either ; a 1D vector or a 2D array of corresponding latitude values. If data is 2D, ; and the LONS parameter is missing, a vector of appropriate length scaled into ; the range -90.0 to 90.0 will be created. ; ; :Keywords: ; cubic: in, optional, type=boolean, default=0 ; If this keyword is set, and the data is a two-dimensional grid, then cubic ; interpolation will be used to create the output image. It is ignored in the ; case of non-gridded input data. ; griddata: in, optional, type=boolean, default=1 ; If the input data is non-gridded, setting this keyword will choose the GRIDDATA ; function to grid the data into a 2D output array. If not set, the data will be ; gridded using the TRIGRID function. ; map: in, optional, type=object ; An input map projection object (cgMap). If provided, the data will be gridded into ; this map projection. If not provided, a map object using a equirectangular map projection ; with a spherical datum will be used. The XRANGE and YRANGE properties of the map object ; will be set by the program in the course of doing the gridding if the `SetRange` keyword is ; set. ; missing: in, optional, type=varies ; Missing data in the gridding process will be set to this value. ; nearest_neighbor: in, optional, type=boolean, default=0 ; If this keyword is set, the nearest neighbor algorithm will be used to create the output ; grid. Otherwise, bilinear (gridded input data) or natural neighbor (non-gridded input data) ; interpolation is used as the default algorithm. The keyword is ignored if non-grided input ; data is being used and the GRIDDATA keyword is not set, or if gridded input data is being ; used and the CUBIC keyword is set. ; resolution: in, optional, type=integer ; A two-element array giving the pixel resolution of the output array in X and Y. ; The default is a 400x400 array. ; setrange: in, optional, type=boolean, default=1 ; If this keyword is set, the XRANGE and YRANGE parameters of the cgMap object will ; be set to the output X and Y ranges. ; xrange: out, optional, type=float ; The output X range in projected meter space (usually associated with the longitude). ; yrange: out, optional, type=float ; The output Y range in projected meter space (usually associated with the latitude). ; ; :Examples: ; To display a GOES image with map annotations:: ; fileURL = 'http://www.idlcoyote.com/misc/goes_example_data.sav' ; filename = "goes_example_data.sav" ; netObject = Obj_New('IDLnetURL') ; void = netObject -> Get(URL=fileURL, FILENAME=filename) ; Obj_Destroy, netObject ; Restore, filename ; peru_lat = Temporary(peru_lat) / 10000. ; peru_lon = Temporary(peru_lon) / 10000. ; s = Size(peruimage, /DIMENSIONS) ; centerLat = peru_lat[s[0]/2, s[1]/2] ; centerLon = peru_lon[s[0]/2, s[1]/2] ; map = Obj_New('cgMap', 'Albers Equal Area', Ellipsoid='sphere', /OnImage, $ ; STANDARD_PAR1=-19, STANDARD_PAR2=20, CENTER_LAT=centerLat, CENTER_LON=centerLon) ; warped = cgWarpToMap(peruImage, peru_lon, peru_lat, MAP=map, MISSING=0, $ ; Resolution=[400, 300], /SetRange) ; cgDisplay, /Free, Title='Warped Image with cgWarpToMap' ; cgImage, warped, Stretch=2, Position=[0,0,1,1] ; map -> Draw ; cgMap_Grid, Map=map, /Label, Color='goldenrod' ; cgMap_Continents, MAP=map, Color='goldenrod' ; cgMap_Continents, MAP=map, Color='goldenrod', /Countries ; ; Additional examples can be found here: http://www.idlcoyote.com/map_tips/warptomap.php. ; ; :Author: ; FANNING SOFTWARE CONSULTING:: ; David W. Fanning ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 ; E-mail: david@idlcoyote.com ; Coyote's Guide to IDL Programming: http://www.idlcoyote.com/ ; ; :History: ; Modification History:: ; Written by David W. Fanning, 12 Sept 2012. ; Modifications to accommodate lat/lon arrays that are one-dimensional to go along ; with 2D data. 13 Sept 2012. DWF. ; ; :Copyright: ; Copyright (c) 2012, Fanning Software Consulting, Inc. ;- FUNCTION cgWarpToMap, data, lons, lats, $ CUBIC=cubic, $ GRIDDATA=griddata, $ MAP=map, $ MISSING=missing, $ NEAREST_NEIGHBOR=nearest_neighbor, $ RESOLUTION=resolution, $ SETRANGE=setrange, $ XRANGE=xrange, $ YRANGE=yrange Compile_Opt idl2 Catch, theError IF theError NE 0 THEN BEGIN Catch, /CANCEL void = cgErrorMsg() RETURN, data ENDIF IF N_Params() EQ 0 THEN BEGIN Print, 'Calling Syntax: warpedImage = cgWarpImage(image, lons, lats, MAP=map, RESOLUTION=resolution)' RETURN, 0 ENDIF ; Handle keywords IF N_Elements(missing) EQ 0 THEN missing = 0B IF Keyword_Set(nearest_neighbor) THEN method = "Nearest Neighbor" ELSE method = "Natural Neighbor" ; Assume we are working with a 2D grid. grid = 1 ; Assume the latitude and longitude arrays are 2D. latlon2d = 1 ; If no map object, use a Equirectangular grid with a spherical datum. IF N_Elements(map) EQ 0 THEN map = Obj_New('cgMap', 'Equirectangular', ELLIPSOID='Sphere') ; Is the data 1D or 2D. If 1D, lons and lats are required. ndim = Size(data, /N_DIMENSIONS) CASE ndim OF 1: BEGIN IF (N_Elements(lons) EQ 0) || (N_Elements(lats) EQ 0) THEN BEGIN Message, 'Must supply longitudes and latitudes with 1D input data.' ENDIF grid = 0 latlon2d = 0 END 2: BEGIN s = Size(data, /DIMENSIONS) IF N_Elements(lons) EQ 0 THEN lons = cgScaleVector(Findgen(s[0]), -180, 180) IF N_Elements(lats) EQ 0 THEN lats = cgScaleVector(Findgen(s[1]), -90, 90) IF Size(lons, /N_DIMENSIONS) EQ 1 THEN latlon2d = 0 IF Size(lats, /N_DIMENSIONS) EQ 1 THEN latlon2d = 0 END ELSE: Message, 'Input data must be either 1D or 2D.' ENDCASE ; Make sure the longitudes are in the range -180 to 180. lons = ((lons + 180) MOD 360) - 180 ; Convert to XY projected meter space. IF N_Elements(lons) NE N_Elements(lats) THEN BEGIN s = Size(data, /DIMENSIONS) lattemp = Rebin(Reform(lats, 1, s[1]), s[0], s[1]) lontemp = Rebin(lons, s[0], s[1]) xy = map -> Forward(lontemp, lattemp) ENDIF ELSE BEGIN xy = map -> Forward(lons, lats) ENDELSE x = Reform(xy[0,*]) y = Reform(xy[1,*]) xmin = Min(x, MAX=xmax) ymin = Min(y, MAX=ymax) ; Set the output resolution of the grid. IF N_Elements(resolution) EQ 0 THEN BEGIN delta_x = (xmax - xmin) / 399 delta_y = (ymax - ymin) / 399 resolution = [400,400] ENDIF ELSE BEGIN delta_x = (xmax - xmin) / (resolution[0]-1) delta_y = (ymax - ymin) / (resolution[1]-1) ENDELSE ; Find the sides of the image and make a boundary rectangle. IF (grid EQ 0) || (latlon2d EQ 0) THEN BEGIN xmin = Min(x, MAX=xmax) ymin = Min(y, MAX=ymax) ENDIF ELSE BEGIN dims = Size(data, /DIMENSIONS) x = Reform(x, dims[0], dims[1]) y = Reform(y, dims[0], dims[1]) xmin = x[0,Round(s[1]/2.)] xmax = x[s[0]-1, Round(s[1]/2.)] ymin = y[Round(s[0]/2.), 0] ymax = y[Round(s[0]/2.), s[1]-1] IF ymin GT ymax THEN BEGIN temp = ymin ymin = ymax ymax = temp ENDIF ENDELSE rect = [xmin-(delta_x/2.), ymin-(delta_y/2.), xmax+(delta_x/2.), ymax+(delta_y/2.)] xrange = rect[[0,2]] yrange = rect[[1,3]] IF Keyword_Set(setrange) THEN map -> SetProperty, XRANGE=xrange, YRANGE=yrange ; If you don't have a grid, then you have to do the gridding the slow way. IF grid EQ 0 THEN BEGIN ; We need a set of Delaunay triangles. QHull, x, y, triangles, /Delaunay ; A choice of GridData or TriGrid for the actual gridding. IF Keyword_Set(griddata) THEN BEGIN warpedImage = GridData(x, y, data, TRIANGLES=triangles, DELTA=[delta_x, delta_y], $ DIMENSION=resolution, START=[Min(x), Min(y)], MISSING=missing, METHOD=method) ENDIF ELSE BEGIN warpedImage = TriGrid(x, y, data, triangles, [delta_x, delta_y], rect, MISSING=missing) ENDELSE RETURN, warpedImage ENDIF ; If you have a grid, you can do the gridding the fast way by interpolating the output ; grid from the input grid. First, create an output grid. xvec = cgScaleVector(Findgen(resolution[0]), xmin-(delta_x/2.), xmax+(delta_x/2.)) yvec = cgScaleVector(Findgen(resolution[1]), ymin-(delta_y/2.), ymax+(delta_x/2.)) x_out = Rebin(xvec, resolution[0], resolution[1]) y_out = Rebin(Reform(yvec, 1, resolution[1]), resolution[0], resolution[1]) ; Get the fractional indices of the output grid on the input grid. dims = Size(data, /DIMENSIONS) xindex = cgScaleVector(x_out, 0, dims[0], Min=xmin, Max=xmax) yindex = cgScaleVector(y_out, 0, dims[1], MIN=ymin, MAX=ymax) ; Interpolate the data. Nearest neighbor, bilinear, or cubic interpolation is possible. IF Keyword_Set(cubic) THEN BEGIN warpedImage = Interpolate(data, xindex, yindex, CUBIC=-0.5, MISSING=missing) ENDIF ELSE BEGIN IF Keyword_Set(nearest_neighbor) THEN BEGIN warpedImage = data[Round(xindex), Round(yindex)] ENDIF ELSE BEGIN warpedImage = Interpolate(data, xindex, yindex, MISSING=missing) ENDELSE ENDELSE RETURN, warpedImage END