File: cggeomosaic.pro

package info (click to toggle)
coyote 2019.01.29-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 6,316 kB
  • sloc: python: 184; makefile: 14; sh: 13
file content (293 lines) | stat: -rw-r--r-- 13,663 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
; docformat = 'rst'
;
; NAME:
;   cgGeoMosaic
;
; PURPOSE:
;   Creates a mosaic or combination image, given the names of two GeoTiff image files.
;   The images must be 2D images.
;
;******************************************************************************************;
;                                                                                          ;
;  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.                            ;
;******************************************************************************************;
;
;+
; Creates a mosaic or combination image, given the names of two GeoTiff image files.
; The images must be 2D images.
; 
; :Categories:
;    Map Utilities
;    
; :Returns:
;    The name of a GeoTiff image file containing the combined image mosaic.
;    
; :Params:
;    geofile_1: in, required, type=string
;        The name of the first GeoTiff file to be combined into a mosaic.
;        Currently the code only works with 2D images.
;    geofile_2: in, required, type=string
;        The name of the second GeoTiff file to be combined into a mosaic.
;        Currently the code only works with 2D images.
;       
; :Keywords:
;     filename: in, optional, type=string
;         The name of the output GeoTiff file. If not provided, the user will
;         be asked to create a name. If provided, the user will not be prompted
;         and this name will be used.
;     imageout: out, optional, type=varies
;         The final image mosaic.
;     mapout: out, optional, type=object
;         A map coordinate object (cgMap) that geolocates the new image mosaic.
;     missing: in, optional, type=varies
;         The missing value in the input images. If scalar value, the same value is used
;         for both images, but may be a two-element array. 
;     miss_out: in, optional, type=varies
;         The missing value of the output image. If undefined, missing[0] is used. If
;         missing[0] is undefined, the value 0B is used.
;     palette: in, optional
;         A 256x3 array containing the R, G, and B color vectors with which the image
;         should be displayed. The palette will be stored with the mosaiced image. If
;         this is not supplied, and either of the other images have a palette, it will be
;         used instead. Images are checked for palettes in the order they are supplied
;         to the function.
;     success: out, optional, type=boolean
;         This keyword is set to 1 if the function completed successfully. And to
;         0 otherwise.
;                
; :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, 18 August 2012. DWF. 
;        Added SUCCESS keyword 4 September 2012. DWF.
;        Now blending overlap regions using 50% of pixel values from the two images. 14 Sept 2012. DWF.
;        Revamp of algorithm's handing of missing values.  Added MISS_OUT keyword and removed HISTOMATCH
;           keyword because it only works for BYTE data. Restricted mosaicking to 2D images. 29 Jan 2013. DWF.
;        Added PALETTE keyword and changed algoritm to use color palettes if available. 21 Nov 2013. DWF.
;
; :Copyright:
;     Copyright (c) 2012-2013, Fanning Software Consulting, Inc.
;-
FUNCTION cgGeoMosaic, geofile_1, geofile_2, $
   FILENAME=filename, $
   IMAGEOUT=newImage, $
   MAPOUT=mapout, $
   MISSING=missing, $
   MISS_OUT=miss_out, $
   PALETTE=palette, $
   SUCCESS=success

    Compile_Opt idl2
    
    ; Error handling
    Catch, theError
    IF theError NE 0 THEN BEGIN
        Catch, /CANCEL
        void = cgErrorMsg()
        IF Obj_Valid(m1) THEN Obj_Destroy, m1
        IF Obj_Valid(m2) THEN Obj_Destroy, m2
        success = 0
        RETURN, ""
    ENDIF
    
    ; Need two files to do the job!
    IF N_Params() NE 2 THEN BEGIN
        void = Dialog_Message('Calling Syntax: geofile = cgGeoMosaic( geofile_1, geofile_2 )')
        RETURN, ""
    ENDIF
    
    ; Assume success.
    success = 1
    
    ; Make sure the missing values are in two-element array, if needed.
    IF N_Elements(missing) NE 0 THEN BEGIN
       IF N_Elements(missing) EQ 1 THEN missing = [missing, missing]
    ENDIF
    IF N_Elements(miss_out) EQ 0 THEN BEGIN
       IF N_Elements(missing) NE 0 THEN miss_out = missing[0] ELSE miss_out = 0B
    ENDIF
    
    ; Let's read the images and get the bounaries of the map projection.
    m1 = cgGeoMap(geofile_1, IMAGE=image_1, PALETTE=pal_1, BOUNDARY=b1, $
        GEOTIFF=geo1, SUCCESS=success)
        
    ; If you haven't been given a palette, then use this palette, if there is one.
    IF N_Elements(palette) EQ 0 THEN IF N_Elements(pal_1) NE 0 THEN palette = pal_1
    IF ~success THEN RETURN, ""
    
    m2 = cgGeoMap(geofile_2, IMAGE=image_2, PALETTE=pal_2, BOUNDARY=b2, $
        GEOTIFF=geo2, SUCCESS=success)
    IF N_Elements(palette) EQ 0 THEN IF N_Elements(pal_2) NE 0 THEN palette = pal_2
    IF ~success THEN RETURN, ""
        
    ; Check to be sure projection and ellipsoid are the same.
    m1 -> GetProperty, Map_Projection=projection_1, Ellipsoid=ellipsoid_1
    m2 -> GetProperty, Map_Projection=projection_2, Ellipsoid=ellipsoid_2
    IF projection_1 NE projection_2 THEN $
       Message, 'The map projections of the two files are not the same.'
    IF ellipsoid_1 NE ellipsoid_2 THEN $
       Message, 'The map projection ellipsoids of the two files are not the same.'
    
    ; If the projected coordinate system is not the same, we have problems.
    IF geo1.PROJECTEDCSTYPEGEOKEY NE geo2.PROJECTEDCSTYPEGEOKEY THEN BEGIN
        Print, 'File 1: ', geofile_1
        Print, 'File 2: ', geofile_2
        Message, 'The projected coordinate systems of the two files are not the same: ' + $
          StrTrim(geo1.PROJECTEDCSTYPEGEOKEY, 2) + ' and ' + StrTrim(geo2.PROJECTEDCSTYPEGEOKEY,2)
    ENDIF
    
    ; Match histograms? Taken out because it only works for Byte images currently.
    ;IF Keyword_Set(histomatch) THEN image_2 = Histomatch(Temporary(image_2), image_1)
         
    ; The scales can be off by a little. We will use the smaller of the two
    ; to make sure we can accommodate both images in the mosaic.
    xscale_1 = geo1.ModelPixelScaleTag[0]
    yscale_1 = geo1.ModelPixelScaleTag[1]
    xscale_2 = geo2.ModelPixelScaleTag[0]
    yscale_2 = geo2.ModelPixelScaleTag[1]
    xscale = xscale_1 < xscale_2
    yscale = yscale_1 < yscale_2
    
    ; Check for missing values.
    IF N_Elements(missing) NE 0 THEN BEGIN
      missingIndices_1 = Where(image_1 EQ missing[0], missingCnt_1)
      IF missingCnt_1 GT 0 THEN image_1[missingIndices_1] = miss_out
      missingIndices_2 = Where(image_2 EQ missing[1], missingCnt_2)
      IF missingCnt_2 GT 0 THEN image_2[missingIndices_2] = miss_out
    ENDIF

    ; Get the image map ranges and create a mosaic image of the proper size.
    m1 -> GetProperty, XRange=xr_1, YRange=yr_1
    m2 -> GetProperty, XRange=xr_2, YRange=yr_2
    xr = [xr_1[0] < xr_2[0], xr_1[1] > xr_2[1]]
    yr = [yr_1[0] < yr_2[0], yr_1[1] > yr_2[1]]
    
    ; How many pixels are in the new image?
    xNumPixels = Ceil((xr[1] - xr[0]) / xscale)
    yNumPixels = Ceil((yr[1] - yr[0]) / yscale)
    
    ; Create a new image of the proper size. 
    imageType = Size(image_1, /Type)
    newImage = Make_Array( xnumPixels, ynumPixels, TYPE=imageType, VALUE=miss_out)
    
    ; Make a copy of this image. This will be the final image.
    finalImage = newImage
    
    ; Create image extent vectors.
    xvec = cgScaleVector(DIndgen(xnumPixels), xr[0], xr[1])
    yvec = cgScaleVector(DIndgen(ynumPixels), yr[0], yr[1])
    
    ; Locate the position of the first image and position it in the final image.
    xloc = 0 > Value_Locate(xvec, [xr_1[0], xr_1[1]])
    yloc = 0 > Value_Locate(yvec, [yr_1[0], yr_1[1]])
    xsize = xloc[1]-xloc[0]+1
    ysize = yloc[1]-yloc[0]+1 
    finalImage[xloc[0]:xloc[1], yloc[0]:yloc[1]] = Congrid(image_1, xsize, ysize)
    
    ; Locate the position of the second image and position it in the new image.
    xloc = 0 > Value_Locate(xvec, [xr_2[0], xr_2[1]])
    yloc = 0 > Value_Locate(yvec, [yr_2[0], yr_2[1]]) 
    xsize = xloc[1]-xloc[0]+1
    ysize = yloc[1]-yloc[0]+1 
    newImage[xloc[0]:xloc[1], yloc[0]:yloc[1]] = Congrid(image_2, xsize, ysize)
    
    ; Where the final image is missing, and the new image is not missing, move
    ; newImage pixels over to final image.
    movingPixels = Where((finalImage EQ miss_out) AND (newimage NE miss_out), movingCount)
    
    ; If there are overlapping pixels, then take the average value.
    overlap =  Where((finalImage NE miss_out) AND (newimage NE miss_out), count)
    finalImage[overlap] = Convert_to_Type(Float(finalImage[overlap])*0.5 + Float(newImage[overlap])*0.5, Size(image_1, /Type))

    ; Move the pixels over.
    IF movingCount GT 0 THEN finalImage[movingPixels] = newImage[movingPixels]
    
    ; Cleanup
    Undefine, newImage
    newImage = Temporary(finalImage)
    
    ; Clean up objects created in the program.
    Obj_Destroy, m1
    Obj_Destroy, m2
    
    ; Need a GeoTiff filename?
    IF N_Elements(filename) EQ 0 THEN BEGIN
        filename = cgPickfile(FILE='mosaic.tif')
        IF filename EQ "" THEN RETURN, ""
    ENDIF    
    
    ; If you have a palette, get the palette RGB vectors.
    IF N_Elements(palette) NE 0 THEN BEGIN
        red = palette[*,0]
        grn = palette[*,1]
        blu = palette[*,2]
    ENDIF
    
    ; Update the GeoTiff structure for the new file and write the file.
    newGeo = geo1
    newGeo.ModelTiePointTag[3:4] = [xr[0], yr[1]]
    newGeo.ModelPixelScaleTag[0] = xscale
    newGeo.ModelPixelScaleTag[1] = yscale
    CASE imageType OF
        1: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, $
              RED=red, GREEN=grn, BLUE=blu
        2: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /SHORT, /SIGNED, $
              RED=red, GREEN=grn, BLUE=blu
        3: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /LONG, /SIGNED, $
              RED=red, GREEN=grn, BLUE=blu
        4: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /FLOAT, $
              RED=red, GREEN=grn, BLUE=blu
        5: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /DOUBLE, $
              RED=red, GREEN=grn, BLUE=blu
        6: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /COMPLEX, $
              RED=red, GREEN=grn, BLUE=blu
        9: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /DCOMPLEX, $
              RED=red, GREEN=grn, BLUE=blu
       12: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /SHORT, $
              RED=red, GREEN=grn, BLUE=blu
       13: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /LONG, $
              RED=red, GREEN=grn, BLUE=blu
       14: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /L64, /SIGNED, $
              RED=red, GREEN=grn, BLUE=blu
       15: Write_Tiff, filename, Reverse(newimage, 2), GEOTIFF=newGeo, /L64, $
              RED=red, GREEN=grn, BLUE=blu
    ENDCASE
    
    ; Create a map object for the image, if one is needed.
    IF Arg_Present(mapOut) THEN mapOut = cgGeoMap(filename)
    
    ; Return the name of the GeoTiff file.
    RETURN, filename
    
END