; 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