File: cgwarptomap.pro

package info (click to toggle)
coyote 2019.01.29-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 6,316 kB
  • sloc: python: 184; makefile: 14; sh: 13
file content (312 lines) | stat: -rw-r--r-- 14,358 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
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
; 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.
;        Fixed a problems in which projected meter X values were sometimes reversed in order. 31 Oct 2014. 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 
        IF xmin GT xmax THEN BEGIN
            temp = xmin
            xmin = xmax
            xmax = 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