forked from idl-coyote/coyote
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cgchangemapprojection.pro
195 lines (177 loc) · 9.63 KB
/
cgchangemapprojection.pro
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
; docformat = 'rst'
;
; NAME:
; cgChangeMapProjection
;
; PURPOSE:
; This function warps a map projected image from one map projection to another, using
; Map_Proj_Image to do the warping.
;
;******************************************************************************************;
; ;
; 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. ;
;******************************************************************************************;
;
;+
; This function warps a map projected image from one map projection to another, using
; Map_Proj_Image to do the warping. Useful in general, it is used specifically to
; warp map projected images into the Equirectangular map projection with a WGS 84 ellipsoid
; that is used by Google Earth. See the program `cgImage2KML` to see how it is used.
;
; :Categories:
; Map Projections, Utilities
;
; :Returns:
; An image warped into the output map projection is returned.
;
; :Params:
; image: in, required
; Any 2D or true-color image with or without an alpha channel. This image will
; be warped to the output map projection (`mapOut`).
; mapin: in, required, type=object
; A map coordinate object (cgMap) describing the map projection and coordinates of the
; input image that is to be warped.
;
; :Keywords:
; boundary: in, optional, type=dblarr
; A four-element array specifying the Cartesian (XY) coordinates of the input
; image range, in the form [xmin, ymin, xmax, ymax]. If this parameter is not
; present, the boundary will be obtained from the `MapIn` map coordinate object.
; On output, this variable will contain the Cartesian boundary of the warped image.
; bilinear: in, optional, type=boolean, default=0
; Set this keyword to warp the image with bilinear interpolation. The default is
; to do nearest neighbor interpolation.
; latlonbox: out, optional, type=fltarr
; A four-element float array containing the boundaries of the warped image in
; the [north, south, east, west] form preferred by Google Earth. Listed as degrees.
; mapout: in, optional, type=object
; A map coordinate object (cgMap) describing the map projection and coordinates of the
; warped image. This image will be warped into this map projection. The [XY]Range of
; this object will be set to the `XYRange` of the output image.
; mask: out, optional
; Set this keyword equal to a named variable that will contain a byte array of the same
; dimensions as the output image, containing a mask of the “good” values. Values in the output
; image that were set to `Missing` (that is, the values were off the map) will have a mask value
; of 0, while all other mask values will be 1.
; missing: in, optional, type=integer, default=0
; Set this keyword equal to an integer value to be used for pixels that fall outside of
; the valid map coordinates. The default value is 0.
; xyrange: out, optional, type=dblarr
; The Cartesian (XY) coordinates associated with the output image. These are the map image
; boundaries of the output image.
;
; :Examples:
; To prepare a GeoTiff file for creating a KML overlay on Google Earth::
; netObject = Obj_New('IDLnetURL')
; url = 'http://www.idlcoyote.com/data/AF03sep15b.n16-VIg.tif'
; returnName = netObject -> Get(URL=url, FILENAME=AF03sep15b.n16-VIg.tif')
; Obj_Destroy, netObject
; map = cgGeoMap('AF03sep15b.n16-VIg.tif')
; googleMap = Obj_New('cgMap', 'Equirectangular', Ellipsoid='WGS 84')
; warpedImage = cgChangeMapProjection(image, map, MAPOUT=googleMap)
;
; :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:
; Change History::
; Written, 30 October 2012 by David W. Fanning.
; Fixed a problem with a TRANSPOSE command for true-color images. Bad logic. 4 Jan 2012. DWF.
; Fixed a problem in which Map_Proj_Image was not handling true-color image warping correctly
; and also a problem with handling missing values passed from cgImage2KML correctly. 20 Feb 2013. DWF.
; Forgot to specify the default value for MISSING. 25 Feb 2014. DWF.
;
; :Copyright:
; Copyright (c) 2012, Fanning Software Consulting, Inc.
;-
FUNCTION cgChangeMapProjection, image, mapIn, $
BOUNDARY=boundary, $
BILINEAR=bilinear, $
LATLONBOX=latlonbox, $
MAPOUT=mapOut, $
MASK=mask, $
MISSING=missing, $
XYRANGE=xyrange
Compile_Opt idl2
; Error handling.
ON_Error, 2
; An image and an input map coordinate object are required.
IF N_Elements(image) EQ 0 THEN Message, 'An input image is required.'
IF N_Elements(mapIn) EQ 0 THEN Message, 'An input map coordinate object is required.'
; Assign default values.
IF N_Elements(boundary) EQ 0 THEN mapIn -> GetProperty, BOUNDARY=boundary
bilinear = Keyword_Set(bilinear)
IF N_Elements(mapOut) EQ 0 THEN BEGIN
mapOut = Obj_New('cgMap', 'Equirectangular', Ellipsoid='WGS 84')
ENDIF
; If the missing value is a color, change it to a color triple for true-color images.
IF N_Elements(missing) NE 0 THEN BEGIN
IF (Size(image, /N_Dimensions) GT 2) THEN BEGIN
IF Size(missing, /TNAME) EQ 'STRING' THEN missing = Reform(cgColor(missing, /Triple))
ENDIF
ENDIF ELSE BEGIN
IF (Size(image, /N_Dimensions) GT 2) THEN missing = BytArr(3) ELSE missing = 0B
ENDELSE
; Gather information about the input image. Do this differently if this is a 2D
; image as opposed to a true-color image.
dims = Image_Dimensions(image, XIndex=xindex, YIndex=yindex, TRUEINDEX=trueindex, $
XSize=xsize, YSize=ysize, ALPHACHANNEL=alphachannel)
IF trueindex EQ -1 THEN BEGIN
warped = Map_Proj_Image(image, boundary, Image_Structure=mapIn->GetMapStruct(), $
Map_Structure=mapOut->GetMapStruct(), UVRANGE=xyrange, Missing=missing[0], MASK=mask)
; Update the range values in the output map coordinate object.
mapOut -> SetProperty, XRANGE=xyrange[[0,2]], YRANGE=xyrange[[1,3]]
; Return information of interest to user.
mapOut -> GetProperty, LATLONBOX=latlonbox, BOUNDARY=boundary
ENDIF ELSE BEGIN
IF trueIndex LT 2 THEN img = Transpose(image, [xindex, yindex, trueindex])
frames = alphachannel ? 4 : 3
warped = Make_Array(xsize, ysize, frames, TYPE=Size(image, /TYPE))
mapInStruct = mapIn->GetMapStruct()
mapOutStruct = mapOut->GetMapStruct()
; Warp the three image planes separately.
FOR j=0,frames-1 DO BEGIN
warped[*,*,j] = Map_Proj_Image(img[*,*,j], boundary, Image_Structure=mapInStruct, $
Map_Structure=mapOutStruct, UVRANGE=xyrange, Missing=missing[j], Mask=mask)
ENDFOR
; Handle alpha channel separately.
IF frames EQ 4 THEN BEGIN
warped[*,*,3] = Map_Proj_Image(img[*,*,3], boundary, Image_Structure=mapInStruct, $
Map_Structure=mapOutStruct, UVRANGE=xyrange)
ENDIF
; Update the range values in the output map coordinate object.
mapOut -> SetProperty, XRANGE=xyrange[[0,2]], YRANGE=xyrange[[1,3]]
; Return information of interest to user.
mapOut -> GetProperty, LATLONBOX=latlonbox, BOUNDARY=boundary
ENDELSE
RETURN, warped
END