; docformat = 'rst' ; ; PURPOSE: ; The purpose of this function is to translate a GEOTIFF structure ; (as returned by QUERY_TIFF or READ_TIFF) into a map coordinate ; object (cgMap) that can be used to georeference images with a ; map data coordinate system. The Coyote Library is required. ; ;******************************************************************************************; ; ; ; Copyright (c) 2011, 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 translate a GEOTIFF structure ; (as returned by QUERY_TIFF or READ_TIFF) into a map coordinate ; object (cgMap) that can be used to georeference images with a ; map data coordinate system. The Coyote Library is required. ; ; It is not possible to have a one-to-one mapping between every GeoTIFF file and a map projection ; in IDL, since IDL has a limited number of map projections and datums available. And, even at that, ; I have not implemented all of IDL's map projections or datums, only those that I thought were most ; likely to be encountered in my own work. If you run into a GeoTIFF file that does not work in this ; code (either because of an error in the code or because it is not supported), please contact me. ; I am interested in supporting as many GeoTIFF files as possible and I will take pains to do so if ; I know they are needed. ; ; :Categories: ; Graphics, Map Projections ; ; :Returns: ; A cgMap coordinate object is returned containing the map projection information ; for the image. ; ; :Params: ; image: in, optional, type=various ; A GeoTIFF image. This can optionally be the name of a GeoTiff file. ; If a filename is used, do not pass the geoTiff parameter, as this ; parameter will be obtained from the GeoTiff file. ; ; :Keywords: ; boundary: out, optional, type=array ; A four-element array giving the boundaries of the map projection in the form ; [x0,y0,x1,y1]. This is a more convenient way of expressing the range ; of the map space. ; ccolor: in, optional, type=string, default='Charcoal' ; The name of a color the map continents should be displayed with. The default ; is "charcoal". Color names are those supported by cgColor. ; clip: in, optional, type=boolean, default=0 ; Set this keyword to display the image (assumes the `Display` keyword is set) ; with a two-precent histogram clipping. ; continents: in, optional, type=boolean, default=0 ; If a cgMap object is made successfully, then setting this keyword ; will add a cgMapContinents object to the cgMap object. ; display: in, optional, type=boolean, default=0 ; Set this keyword to display the image with the map annotations in a ; resizeable cgWindow. ; ellipsoid: out, optional, type=string ; The name of the ellipsoid used in the GeoTiff file. ; gcolor: in, optional, type=string, default='Gray' ; The name of a color the map grid should be displayed with. The default ; is "gray". Color names are those supported by cgColor. ; geotiff: in, optional, type=structure ; A GeoTIFF structure of geoTags. Normally obtained by calling ; QUERY_TIFF or READ_TIFF, but will be read from the file if a ; filename is passed as the first positional parameter. If the geotag ; is read from the file, it can be passed back to the user by setting ; this keyword to a named variable. ; grid: in, optional, type=boolean, default=0 ; If a cgMap object is made successfully, then setting this keyword ; will add a cgMapGrid object to the cgMap object. ; image: out, optional, type=varies ; Set this keyword to a named variable that on exit will contain the image data. ; latlonbox: out, optional, type=array ; A four-element array giving the boundaries of the map projection in the ; Google Map form of [north, south, east, west]. This is useful when you ; are creating image overlays to be added to Google Earth. ; map_projection: out, optional, type=string ; The name of the map projection found in the GeoTiff file. ; mcolor: in, optional, type=string, default='Black' ; The name of a color the map should be displayed in. (Normally the map ; border and map title are displayed in this color.) ; onimage: in, optional, type=boolean, default=0 ; Set this keword if the map object is to get its position from the last ; cgImage command issued. ; palette: out, optional, type=bytarr ; If the GeoTiff file contains RGB color vectors, and keywords to cgGeoTiff cause the ; file to be read (e.g, IMAGE or DISPLAY), then this output keyword will contain those ; vectors in a 3-by-256 byte array. ; silent: in, optional, type=boolean, default=0 ; IDL cannot map every GeoTiff image to a supported map projection or datum. ; Normally, if the GeoTIFF image is unsupported, an error message is issued. ; Setting this keyword will suppress such error messages. If you do this, you ; MUST check the SUCCESS keyword to see if the program ran successfully. (Of ; course, it is a good idea to check it in any case!) ; sub_rect: in, optional, type=integer ; Set this keyword to a four-element array, [x, y, width, height], that ; specifies a rectangular region within the file to extract. Only the ; rectangular portion of the image selected by this keyword is read and ; returned. The rectangle is measured in pixels from the lower left ; corner (right hand coordinate system). If this keyword is not use, the ; entire image is read. ; success: out, optional, type=boolean, default=0 ; An output variable that will contain a 1 if the map coordinate object was ; created successfully. Or to a 0 if it was not created successfully. ; title: in, optional, type=string, default="" ; The title of the map projection. ; _extra: in, optional ; Any keyword appropriate for cgImage can be collected and passed along if ; the DISPLAY keyword is also set. ; ; :Examples: ; To display a GeoTiff image in cgWindow with map annotations:: ; file = 'C:\IDL\data\tiff\AF03sep15b.n16-VIg.tif' ; mapCoord = cgGeoMap(file, /Continents, /Grid, /Display) ; ; :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:: ; Converted from old GeoTIff program in Catalyst Library to run with Coyote Graphics routines. ; 9 November 2011. David W. Fanning. ; Assuming meters as unit length if can't find linear measure in file. 1 Dec 2011. DWF. ; Added UTM Zone 18S (projection 29178) for LandSAT MSS images over Amazon. 1 Dec 2011. DWF. ; Added the ability to pass cgImage keywords via the keyword inheritance mechanism. 1 Dec 2011. DWF. ; Added SUB_RECT keyword. 5 Dec 2011. DWF. ; Added ability to read GOES EAST and GOES WEST GeoTiff files from this NOAA web site: ; http://http://www.osdpd.noaa.gov/ml/gis/index.html. 26 Dec 2011. DWF. ; Modified code to read GeoTiff files created by HEG v2.11 from HDF-EOS2 grid files. 30 Dec 2011. DWF. ; Modified to read multi-dimensional GeoTiff images and reverse the images correct. 12 Jan 2012. DWF. ; Write base filename as title for window if the DISPLAY keyword is set. 22 Feb 2012. DWF. ; Had inexplicably left out CENTER_LATITUDE parameter in Equirectangular projection. 30 July 2012. DWF. ; Added PALETTE keyword to return the RGB color palette present in the file, if any. 16 August 2012. DWF. ; Moved the geotiff argument to a GEOTIFF keyword, as I always expect it to be this way. 6 Sept 2012. DWF. ; Removed UTM/WGS84 warning message in IDL 8.2, as this problem has been fixed in IDL 8.2. 2 Oct 2012. DWF. ; Added BOUNDARY, ELLIPSOID, LATLONBOX, and MAP_PROJECTION output keywords to facilitate ; creating Google Earth overlays. 30 Oct 2012. DWF. ; Small fix to allow true-color TIFF images to be read correctly. 11 Sept 2013. DWF. ; Modified program to look for NSIDC Hughes ellipsoid specified as semi-major axis of 6378273.0 and ; semi-minor axis of 6356889.4 as specified by the GeoTiff tags GEOGSEMIMAJORAXISGEOKEY and ; GEOGSEMIMINORAXISGEOKEY. I do this ONLY for Polar Stereographic map projections. ; True-color geoTiff files not being returned correctly because of a typo. Changed "GT" to "GE" on line 1380. 8 Dec 2013. DWF. ; Ran into GeoTiff file with a ModelTransformationTag field instead of the usual ModelPixelScaleTag ; and ModelTiePointTag fields. Modified the program to get needed information from the ; transformation matrix. This *assumes* no rotation, etc. 10 Oct 2014. DWF. ; ; :Copyright: ; Copyright (c) 2011-2014, Fanning Software Consulting, Inc. ;- Function cgGeoMap, image, $ BOUNDARY=boundary, $ CLIP=clip, $ DISPLAY=display, $ ELLIPSOID=ellipsoid, $ GCOLOR=gcolor, $ GEOTIFF=geotiff, $ GRID=grid, $ CCOLOR=ccolor, $ CONTINENTS=continents, $ IMAGE=outImage, $ LATLONBOX=latlonbox, $ MAP_PROJECTION=map_projection, $ MCOLOR=mcolor, $ ONIMAGE=onimage, $ PALETTE=palette, $ SILENT=silent, $ SUB_RECT=sub_rect, $ SUCCESS=success, $ TITLE=title, $ _EXTRA=extra ; Compiler options. Compile_Opt DEFINT32 Compile_Opt STRICTARR Compile_Opt STRICTARRSUBS Compile_Opt LOGICAL_PREDICATE Catch, theError IF theError NE 0 THEN BEGIN Catch, /CANCEL IF ~Keyword_Set(silent) THEN void = cgErrorMsg() success = 0 RETURN, Obj_New() ENDIF ; Set default values of keywords. SetDefaultValue, mcolor, 'black' SetDefaultValue, gcolor, 'charcoal' SetDefaultValue, ccolor, 'charcoal' SetDefaultValue, title, "" onimage = Keyword_Set(onimage) ; If no parameters, ask about reading a GeoTiff file. IF N_Elements(image) EQ 0 THEN BEGIN geofile = cgPickfile(Title='Select GeoTiff File...', FILTER=['*.tif', '*.tiff']) IF geofile EQ "" THEN RETURN, Obj_New() ; Read the image query file. Make sure it is a GeoTiff file. ok = Query_Tiff(geofile, info, GEOTIFF=geotiff) IF ~ok THEN Message, 'Selected file does not appear to be a TIFF file.' IF Size(geotiff, /TNAME) NE 'STRUCT' THEN Message, 'Selected file does not appear to be a GeoTiff file.' channels = info.channels dims = info.dimensions CASE channels OF 1: BEGIN xsize = dims[0] ysize = dims[1] END 3: BEGIN indices = Where(dims NE 3, count) IF count NE 2 THEN Message, 'Cannot determine the image dimensions.' xsize = dims[indices[0]] ysize = dims[indices[1]] END ELSE: Message, 'Do not know how to handle an image with ' + StrTrim(channels,2) + ' channels' ENDCASE ; For further processing we need image to be equal to geofile. image = geofile ENDIF ; Is the first parameter the name of a file or an image? IF (N_Elements(image) NE 0) || (N_Elements(geofile) NE 0) THEN BEGIN ; Is the parameter a string? IF Size(image, /TNAME) NE 'STRING' THEN BEGIN ; Then we are going to assume this is an image and that you have passed ; a geotiff structure to us. IF N_Elements(geotiff) EQ 0 THEN Message, 'A GeoTiff structure must be passed with the image' ; Get the image sizes dims = Image_Dimensions(image, XSIZE=xsize, YSIZE=ysize) ENDIF ELSE BEGIN ; Is it a filename? Can we find the file? ok = File_Test(image) IF ~ok THEN Message, 'Cannot locate the specified file: ' + image ; Read the image query file. Make sure it is a GeoTiff file. ok = Query_Tiff(image, info, GEOTIFF=geotiff) IF ~ok THEN Message, 'Specified file does not appear to be a TIFF file.' IF Size(geotiff, /TNAME) NE 'STRUCT' THEN Message, 'Specified file does not appear to be a GeoTiff file.' channels = info.channels dims = info.dimensions CASE channels OF 1: BEGIN xsize = dims[0] ysize = dims[1] END 3: BEGIN indices = Where(dims NE 3, count) IF count NE 2 THEN Message, 'Cannot determine the image dimensions.' xsize = dims[indices[0]] ysize = dims[indices[1]] END 4: BEGIN indices = Where(dims NE 4, count) IF count NE 2 THEN Message, 'Cannot determine the image dimensions.' xsize = dims[indices[0]] ysize = dims[indices[1]] END ELSE: Message, 'Do not know how to handle an image with ' + StrTrim(channels,2) + ' channels' ENDCASE ENDELSE ENDIF ; Get the fields of the geotiff structure. fields = Tag_Names(geotiff) ; We can only handle raster images with projected coordinate systems, unless this is ; a GeoTiff file with Geographic model, such as GOES images. If it is one of these images, ; then I am going to assign a project code to this image and handle it in the special projection ; section below. geoModel = 0 gtModelIndex = Where(fields EQ 'GTMODELTYPEGEOKEY', gtModelType) IF gtModelType GT 0 THEN BEGIN IF (geotiff.gtModelTypeGeoKey EQ 2) THEN BEGIN ; We need a map projection code if we have one in the file index = Where(fields EQ 'PROJECTEDCSTYPEGEOKEY', projCount) IF projCount GT 0 THEN projCode = geotiff.PROJECTEDCSTYPEGEOKEY index = Where(fields EQ 'PROJECTIONGEOKEY', projCount) IF projCount GT 0 THEN projCode = geotiff.PROJECTIONGEOKEY IF N_Elements(projCode) EQ 0 THEN projCode = 99901 ; GOES East/West Satellites geoModel = 1 ENDIF ENDIF ; We want all measurements to be in meters. First, check to see of there is a ProjLineearUnitsGeokey. ; If not, check the GeogLinearUnitsGeokey. index = Where(fields EQ 'PROJLINEARUNITSGEOKEY', unitsCount) IF unitsCount EQ 0 THEN BEGIN index = Where(fields EQ 'GEOGLINEARUNITSGEOKEY', unitsCount) IF unitsCount NE 0 THEN unitValue = geotiff.GEOGLINEARUNITSGEOKEY ENDIF ELSE unitValue = geotiff.PROJLINEARUNITSGEOKEY ; Can't find a linear unit. Will assume meters. No trouble with this so far. (Fingers crossed.) IF N_Elements(unitValue) EQ 0 THEN BEGIN unitValue = 9001L ; IF ~Keyword_Set(silent) THEN $ ; Message, 'Cannot find a linear unit key in the GeoTiff structure. Assuming meter.', /Informational ENDIF CASE unitValue OF '9001': ; Linear meter '9002': BEGIN ; Linear foot xscale = xscale * (1.0D/3.28D) yscale = yscale * (1.0D/3.28D) END ELSE: Message, 'Linear unit ' + StrTrim(unitValue,2) + ' is not supported in this program.' ENDCASE ; We need a map projection code if we have one in the file index = Where(fields EQ 'PROJECTEDCSTYPEGEOKEY', projCount) IF projCount EQ 0 THEN BEGIN index = Where(fields EQ 'PROJECTIONGEOKEY', projCount) IF projCount GT 0 THEN projCode = geotiff.PROJECTIONGEOKEY ENDIF ELSE projCode = geotiff.PROJECTEDCSTYPEGEOKEY IF N_Elements(projCode) EQ 0 THEN Message, 'Cannot find a map projection Geokey in the GeoTiff Structure.' ; Assume successful completion. success = 1 ; Need to find the pixel scales and the location of the (0,0) point in the ; image. Some images store this information is a MODELTRANSFORMATIONTAG field, ; and others store it in ModelPixelScaleTag and ModelTiePointTag fields. These ; two methods are mutually exclusive. index = Where(fields EQ 'MODELTRANSFORMATIONTAG', modelCount) IF modelCount GT 0 THEN BEGIN ; Get the pixel scale values. xscale = (geotiff.MODELTRANSFORMATIONTAG)[0,0] yscale = Abs((geotiff.MODELTRANSFORMATIONTAG)[1,1]) ; Get the tie points and calculate the map projection range. x0 = (geotiff.MODELTRANSFORMATIONTAG)[3,0] y1 = (geotiff.MODELTRANSFORMATIONTAG)[3,1] ENDIF ELSE BEGIN ; Get the pixel scale values. xscale = (geotiff.ModelPixelScaleTag)[0] yscale = (geotiff.ModelPixelScaleTag)[1] ; Get the tie points and calculate the map projection range. x0 = (geotiff.ModelTiePointTag)[3] y1 = (geotiff.ModelTiePointTag)[4] ENDELSE x1 = x0 + (xscale * xsize) y0 = y1 - (yscale * ysize) xrange = [x0, x1] yrange = [y0, y1] ; Print, 'X Range: ', xrange ; Print, 'Y Range: ', yrange ; Is this a user-defined projection code? IF projCode EQ 32767L THEN BEGIN ; What kind of datum is being used? IF Where(fields EQ 'GEOGRAPHICTYPEGEOKEY') NE -1 THEN BEGIN CASE geotiff.GEOGRAPHICTYPEGEOKEY OF 0: thisDatum = 19 ; Undefined, assgined to sphere. 4001: thisDatum = 9 ; Airy 4002: thisDatum = 11 ; Modified Airy 4003: thisDatum = 14 ; Astralian National 4008: thisDatum = 0 ; Clarke 1866 4017: thisDatum = 6 ; Everest 4018: thisDatum = 10 ; Modified Everest 4023: thisDatum = 3 ; International 1967 4024: thisDatum = 15 ; Krassovsky 4030: thisDatum = 8 ; WGS 84 4034: thisDatum = 1 ; Clarke 1880 4035: thisDatum = 19 ; Sphere 4054: thisDatum = 25 ; Hughes 4267: thisDatum = 0 ; NAD27, same as Clark 1866 4322: thisDatum = 5 ; WGS 72 4326: thisDatum = 8 ; WGS 84 ELSE: ENDCASE ENDIF ; What kind of datum ellipsoid is being used? IF Where(fields EQ 'GEOGELLIPSOIDGEOKEY') NE -1 THEN BEGIN CASE geotiff.GEOGELLIPSOIDGEOKEY OF 0: thisDatum = 19 ; Undefined, assgined to sphere. 6322: thisDatum = 5 ; WGS 72 6326: thisDatum = 8 ; WGS 84 6001: thisDatum = 9 ; Airy 6002: thisDatum = 11 ; Modified Airy 6003: thisDatum = 14 ; Astralian National 6008: thisDatum = 0 ; Clarke 1866 6017: thisDatum = 6 ; Everest 6018: thisDatum = 10 ; Modified Everest 6023: thisDatum = 3 ; International 1967 6024: thisDatum = 15 ; Krassovsky 6030: thisDatum = 8 ; WGS 84 6034: thisDatum = 1 ; Clarke 1880 6035: thisDatum = 19 ; Sphere ELSE: ENDCASE ENDIF ; What kind of projection is this? IF Where(fields EQ 'PROJCOORDTRANSGEOKEY') NE -1 THEN BEGIN CASE geotiff.PROJCOORDTRANSGEOKEY OF 1: thisProjection = 109 ; Transverse Mercator 7: thisProjection = 105 ; Mercator 8: thisProjection = 104 ; Lambert Conformal Conic 10: thisProjection = 111 ; Lambert Azimuthal Equal Area 11: thisProjection = 103 ; Albers Equal Area 14: thisProjection = 110 ; Stereographic 15: thisProjection = 106 ; Polar Stereographic 17: thisProjection = 117 ; Equirectangular 20: thisProjection = 118 ; Miller Cylindrical 21: thisProjection = 114 ; Orthographic 22: thisProjection = 107 ; Polyconic 23: thisProjection = 121 ; Robinson 24: thisProjection = 116 ; Sinusoidal ELSE: Message, 'Unsupported map projection for this GeoTIFF image.' ENDCASE ENDIF IF Where(fields EQ 'PROJCOORDTRANSGEOKEY') NE -1 THEN BEGIN CASE geotiff.PROJCOORDTRANSGEOKEY OF 1: thisProjection = 109 ; Transverse Mercator 7: thisProjection = 105 ; Mercator 8: thisProjection = 104 ; Lambert Conformal Conic 10: thisProjection = 111 ; Lambert Azimuthal Equal Area 11: thisProjection = 103 ; Albers Equal Area 14: thisProjection = 110 ; Stereographic 15: thisProjection = 106 ; Polar Stereographic 17: thisProjection = 117 ; Equirectangular 20: thisProjection = 118 ; Miller Cylindrical 21: thisProjection = 114 ; Orthographic 22: thisProjection = 107 ; Polyconic 23: thisProjection = 121 ; Robinson 24: thisProjection = 116 ; Sinusoidal ELSE: Message, 'Unsupported map projection for this GeoTIFF image.' ENDCASE ENDIF ELSE BEGIN IF geoModel THEN thisProjection = 117 ENDELSE IF N_Elements(thisProjection) EQ 0 THEN Message, 'Unknown map projection for this GeoTIFF image.' ; To accommodate NSIDC polar stereographic projection with a HUGHES ellipsoid, look for length of semi-major and ; semi-minor axes if datum is undefined and this is a Polar Stereographic map projection. IF (thisProjection EQ 106) && (N_Elements(thisDatum) EQ 0) THEN BEGIN IF (Where(fields EQ 'GEOGSEMIMAJORAXISGEOKEY'))[0] NE -1 THEN semimajor = geotiff.GEOGSEMIMAJORAXISGEOKEY IF (Where(fields EQ 'GEOGSEMIMINORAXISGEOKEY'))[0] NE -1 THEN semiminor = geotiff.GEOGSEMIMINORAXISGEOKEY IF (N_Elements(semimajor) NE 0) AND (N_Elements(semiminor) NE 0) THEN BEGIN IF (semimajor GE 6378272.0d) && (semimajor LE 6378274.0d) && $ (semiminor GE 6356889.3d) && (semiminor LE 6356889.5d) THEN thisDatum = 26 ; Hughes ellipsoid ENDIF ENDIF ; Default datum if you don't have one defined by here. IF N_Elements(thisDatum) EQ 0 THEN thisDatum = 19 ; Defaults to Sphere. ; Make the map coordinate object, based on the map projection. CASE thisProjection OF 103: BEGIN ; Albers Equal Area CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 std_parallel_1 = geotiff.PROJSTDPARALLEL1GEOKEY std_parallel_2 = geotiff.PROJSTDPARALLEL2GEOKEY mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ STANDARD_PAR1=std_parallel_1, STANDARD_PAR2=std_parallel_2, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, NAME='FROM_GEOCOORD', $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 104: BEGIN ; Lambert Conformal Conic CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 std_parallel_1 = geotiff.PROJSTDPARALLEL1GEOKEY std_parallel_2 = geotiff.PROJSTDPARALLEL2GEOKEY mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ STANDARD_PAR1=std_parallel_1, STANDARD_PAR2=std_parallel_2, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 105: BEGIN ; Mercator CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 index = Where(fields EQ 'PROJSTDPARALLEL1GEOKEY', count) IF count EQ 0 THEN BEGIN index = Where(fields EQ 'PROJSTDPARALLELGEOKEY', count) IF count NE 0 THEN true_scale_lat = geotiff.PROJSTDPARALLELGEOKEY ENDIF ELSE true_scale_lat = geotiff. PROJSTDPARALLEL1GEOKEY mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ TRUE_SCALE_LATITUDE=true_scale_lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 106: BEGIN ; Polar Stereographic CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 107: BEGIN ; Polyconic CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 109: BEGIN ; Transverse Mercator CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 index = Where(fields EQ 'PROJSCALEATNATORIGINGEOKEY', count) IF count GT 0 THEN mercator_scale = geotiff.PROJSCALEATNATORIGINGEOKEY ELSE mercator_scale = 0.9996 mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ MERCATOR_SCALE=mercator_scale, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 110: BEGIN ; Stereographic CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 111: BEGIN ; Lambert Azimuthal CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 114: BEGIN ; Orthographic CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 116: BEGIN ; Sinusoidal CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 117: BEGIN ; Equirectangular CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 index = Where(fields EQ 'PROJSTDPARALLEL1GEOKEY', count) IF count EQ 0 THEN BEGIN index = Where(fields EQ 'PROJSTDPARALLELGEOKEY', count) IF count NE 0 THEN true_scale_lat = geotiff.PROJSTDPARALLELGEOKEY ENDIF ELSE true_scale_lat = geotiff. PROJSTDPARALLEL1GEOKEY mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ TRUE_SCALE_LATITUDE=true_scale_lat, CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 118: BEGIN ; Miller Cylindrical CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END 121: BEGIN ; Robinson CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE index = Where(fields EQ 'PROJFALSEEASTINGGEOKEY', count) IF count GT 0 THEN falseEast = geotiff.PROJFALSEEASTINGGEOKEY ELSE falseEast = 0 index = Where(fields EQ 'PROJFALSENORTHINGGEOKEY', count) IF count GT 0 THEN falseNorth = geotiff.PROJFALSENORTHINGGEOKEY ELSE falseNorth = 0 mapCoord = Obj_New('cgMap', thisProjection, DATUM=19, $ CENTER_LONGITUDE=lon, $ FALSE_EASTING=falseEast, FALSE_NORTHING=falseNorth, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange END ENDCASE ; RETURN, mapCoord ENDIF ; UTM projections. IF (projCode GT 32200) AND (projCode LE 32260) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 5 ; WGS 72 zone = projCode - 32200 CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF IF (projCode GT 32300) AND (projCode LE 32360) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 5 ; WGS 72 zone = -(projCode - 32300) CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF IF (projCode GT 32400) AND (projCode LE 32460) THEN BEGIN Message, 'The WGS72BE datum is not currently supported. ; thisProjection = 101 ; UTM ; thisDatum = 5 ; WGS 72 ; zone = -(projCode - 32400) ; index = Where(fields EQ 'PROJNATORIGINLONGGEOKEY', count) ; IF count EQ 0 THEN BEGIN ; index = Where(fields EQ 'PROJORIGINLONGGEOKEY', count) ; IF count NE 0 THEN lon = geotiff.PROJORIGINLONGGEOKEY ; ENDIF ELSE lon = geotiff.PROJNATORIGINLONGGEOKEY ; index = Where(fields EQ 'PROJNATORIGINLATGEOKEY', count) ; IF count EQ 0 THEN BEGIN ; index = Where(fields EQ 'PROJORIGINLATGEOKEY', count) ; IF count NE 0 THEN lat = geotiff.PROJORIGINLATGEOKEY ; ENDIF ELSE lat = geotiff. PROJNATORIGINLATGEOKEY ; mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ ; CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone) ; mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF IF (projCode GT 32500) AND (projCode LE 32560) THEN BEGIN Message, 'The WGS72BE datum is not currently supported. ; thisProjection = 101 ; UTM ; thisDatum = 5 ; WGS 72 ; zone = -(projCode - 32500) ; index = Where(fields EQ 'PROJNATORIGINLONGGEOKEY', count) ; IF count EQ 0 THEN BEGIN ; index = Where(fields EQ 'PROJORIGINLONGGEOKEY', count) ; IF count NE 0 THEN lon = geotiff.PROJORIGINLONGGEOKEY ; ENDIF ELSE lon = geotiff.PROJNATORIGINLONGGEOKEY ; index = Where(fields EQ 'PROJNATORIGINLATGEOKEY', count) ; IF count EQ 0 THEN BEGIN ; index = Where(fields EQ 'PROJORIGINLATGEOKEY', count) ; IF count NE 0 THEN lat = geotiff.PROJORIGINLATGEOKEY ; ENDIF ELSE lat = geotiff. PROJNATORIGINLATGEOKEY ; mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ ; CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone) ; mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF IF (projCode GT 32600) AND (projCode LE 32660) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 8 ; WGS 84 zone = projCode - 32600 ; There is a bug in all versions of IDL up to IDL 8.1 apparently that ; produces the wrong result when a UTM projection is used in conjunction ; with a WGS84 datum (the most common datum used in this projection). Here ; we substitute the WALBECK datum, which is nearly identical to WGS84 are ; results in position errors of less than a meter typically. IF (Float(!version.release) LT 8.2) THEN BEGIN IF ~Keyword_Set(silent) THEN Print, 'Switching UTM datum from WGS84 to WALBECK to avoid UTM projection bug.' thisDatum = 12 ENDIF CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF IF (projCode GT 32700) AND (projCode LE 32760) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 8 ; WGS 84 zone = -(projCode - 32700) ; There is a bug in all versions of IDL up to IDL 8.1 apparently that ; produces the wrong result when a UTM projection is used in conjunction ; with a WGS84 datum (the most common datum used in this projection). Here ; we substitute the WALBECK datum, which is nearly identical to WGS84 are ; results in position errors of less than a meter typically. IF (Float(!version.release) LT 8.2) THEN BEGIN IF ~Keyword_Set(silent) THEN Print, 'Switching UTM datum from WGS84 to WALBECK to avoid UTM projection bug.' thisDatum = 12 ENDIF CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF ; State Zone Projections IF (projCode GT 26700) AND (projCode LE 26899) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 0 ; CLARK 1866 or NAD27 zone = projCode - 26700 CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF IF (projCode GT 26900) AND (projCode LE 27999) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 8 ; WGS84 or NAD83 zone = projCode - 26900 ; There is a bug in all versions of IDL up to IDL 8.1 apparently that ; produces the wrong result when a UTM projection is used in conjunction ; with a WGS84 datum (the most common datum used in this projection). Here ; we substitute the WALBECK datum, which is nearly identical to WGS84 are ; results in position errors of less than a meter typically. IF (Float(!version.release) LT 8.2) THEN BEGIN IF ~Keyword_Set(silent) THEN Print, 'Switching UTM datum from WGS84 to WALBECK to avoid UTM projection bug.' thisDatum = 12 ENDIF CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ; RETURN, mapCoord ENDIF ; State Zone Projections IF (projCode GT 32000) AND (projCode LE 32099) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 0 ; CLARK 1866 or NAD27 zone = projCode - 32000 CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ENDIF IF (projCode GT 32100) AND (projCode LE 32299) THEN BEGIN thisProjection = 101 ; UTM thisDatum = 8 ; WGS84 or NAD83 zone = projCode - 32100 ; There is a bug in all versions of IDL up to IDL 8.1 apparently that ; produces the wrong result when a UTM projection is used in conjunction ; with a WGS84 datum (the most common datum used in this projection). Here ; we substitute the WALBECK datum, which is nearly identical to WGS84 are ; results in position errors of less than a meter typically. IF (Float(!version.release) LT 8.2) THEN BEGIN IF ~Keyword_Set(silent) THEN Print, 'Switching UTM datum from WGS84 to WALBECK to avoid UTM projection bug.' thisDatum = 12 ENDIF CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ENDIF ; Projections that have been added to support work I am doing for clients. ; LandSAT MSS Amazon Basin IF (projCode EQ 29178) THEN BEGIN ; UTM Zone 18S thisProjection = 101 ; UTM IF (Float(!Version.Release) GE 8) $ THEN thisDatum = 23 $ ; South American 1969 ELSE thisDatum = 14 ; Australian National (same as SAD-69) zone = -18 ; Zone 18S CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ENDIF IF (projCode EQ 29181) THEN BEGIN ; UTM Zone 21S thisProjection = 101 ; UTM IF (Float(!Version.Release) GE 8) $ THEN thisDatum = 23 $ ; South American 1969 ELSE thisDatum = 14 ; Australian National (same as SAD-69) zone = -21 ; Zone 21S CASE 1 OF ; Longitude (Where(fields EQ 'PROJNATORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJNATORIGINLONGGEOKEY (Where(fields EQ 'PROJORIGINLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJORIGINLONGGEOKEY (Where(fields EQ 'PROJCENTERLONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJCENTERLONGGEOKEY (Where(fields EQ 'PROJSTRAIGHTVERTPOLELONGGEOKEY'))[0] NE -1 : lon = geotiff.PROJSTRAIGHTVERTPOLELONGGEOKEY ELSE: ENDCASE CASE 1 OF ; Latitude (Where(fields EQ 'PROJNATORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJNATORIGINLATGEOKEY (Where(fields EQ 'PROJORIGINLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJORIGINLATGEOKEY (Where(fields EQ 'PROJCENTERLATGEOKEY'))[0] NE -1 : lat = geotiff.PROJCENTERLATGEOKEY ELSE: ENDCASE mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LATITUDE=lat, CENTER_LONGITUDE=lon, ZONE=zone, $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage) mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange ENDIF IF (projCode EQ 99901) THEN BEGIN ; GOES EAST image thisProjection = 117 ; Equirectangular ; What kind of datum is being used? IF Where(fields EQ 'GEOGRAPHICTYPEGEOKEY') NE -1 THEN BEGIN CASE geotiff.GEOGRAPHICTYPEGEOKEY OF 0: thisDatum = 19 ; Undefined, assgined to sphere. 4001: thisDatum = 9 ; Airy 4002: thisDatum = 11 ; Modified Airy 4003: thisDatum = 14 ; Astralian National 4008: thisDatum = 0 ; Clarke 1866 4017: thisDatum = 6 ; Everest 4018: thisDatum = 10 ; Modified Everest 4023: thisDatum = 3 ; International 1967 4024: thisDatum = 15 ; Krassovsky 4030: thisDatum = 8 ; WGS 84 4034: thisDatum = 1 ; Clarke 1880 4035: thisDatum = 19 ; Sphere 4267: thisDatum = 0 ; NAD27, same as Clark 1866 4322: thisDatum = 5 ; WGS 72 4326: thisDatum = 8 ; WGS 84 ELSE: ENDCASE ENDIF ; Get the pixel scale values, which are in lat/lon degrees in this case. xscale = (geotiff.ModelPixelScaleTag)[0] yscale = (geotiff.ModelPixelScaleTag)[1] ; Get the tie points (in lat/lon space) and calculate the map projection range. x0 = (geotiff.ModelTiePointTag)[3] y0 = (geotiff.ModelTiePointTag)[4] x1 = x0 + (xscale * xsize) y1 = y0 - (yscale * ysize) xrange = [x0, x1] yrange = [y0, y1] mapCoord = Obj_New('cgMap', thisProjection, DATUM=thisDatum, $ CENTER_LONGITUDE=(xrange[1]-xrange[0])/2.0 + xrange[0], $ CENTER_LATITUDE=(yrange[1]-yrange[0])/2.0 + yrange[0], $ COLOR=mcolor, TITLE=title, ONIMAGE=onimage, XRANGE=xrange, YRANGE=yrange) ENDIF ; If this is a Geographic Model, then the ranges are in lon/lat space. ; They need to be converted to XY projected space. IF geoModel THEN BEGIN mapCoord -> GetProperty, XRANGE=xrange, YRANGE=yrange ; Many of these geographic models go between 0 and 360 degrees ; of longitude, which cause me a great deal of headache. If the ; range is 360 degrees, more or less, then I'm just going to ; assign the range. IF Floats_Equal(xrange[1] - xrange[0], 360.0, ULP=5) THEN xrange = [-180, 180] ; Otherwise, I will project them into XY meters. xy = mapCoord -> Forward(xrange, yrange) mapCoord -> SetProperty, XRANGE=Reform(xy[0,*]), YRANGE=Reform(xy[1,*]) ENDIF ; If we get here, we don't know what to do. IF (N_Elements(mapCoord) EQ 0) THEN BEGIN Message, 'Map projection code ' + StrTrim(projCode,2) + ' is not supported. ENDIF ; Need a map outline in this cgMapCoord object? IF Keyword_Set(continents) THEN BEGIN contObj = Obj_New('cgMapContinents', mapCoord, /HIRES, $ COLOR=ccolor, FILL=Keyword_Set(fill), NAME='MAPCONTINENTS') IF ~Obj_Valid(contObj) THEN $ Message, 'Cannot successfully create a cgMapContinents object.' mapCoord -> SetProperty, CONTINENTS=contObj ENDIF ; Need a map grid in this cgMapCoord object? IF Keyword_Set(grid) THEN BEGIN gridObj = Obj_New('cgMapGrid', mapCoord, COLOR=gcolor, /AUTODRAWGRID, NAME='MAPGRID') IF ~Obj_Valid(gridObj) THEN $ Message, 'Cannot successfully create a cgMapGrid object.' mapCoord -> SetProperty, GRID=gridObj ENDIF ; Would you like to display the image? IF Keyword_Set(display) THEN BEGIN IF N_Elements(geofile) EQ 0 THEN geofile = image thisImage = Read_Tiff(geofile, r, g, b, SUB_RECT=sub_rect) dims = Image_Dimensions(thisImage, YINDEX=yindex) thisImage = Reverse(Temporary(thisImage), yindex+1) IF N_Elements(r) NE 0 THEN palette=[[r],[g],[b]] mapCoord -> SetProperty, ONIMAGE=1 cgWindow, WASPECT=Float(ysize)/xsize, WTitle=cgRootName(geofile) IF Keyword_Set(clip) THEN BEGIN cgImage, thisImage, /Keep_Aspect, /AddCmd, Margin=0.05, $ STRETCH=2, PALETTE=palette, _EXTRA=extra ENDIF ELSE BEGIN cgImage, thisImage, /Keep_Aspect, /AddCmd, Margin=0.05, PALETTE=palette, _EXTRA=extra ENDELSE cgWindow, 'Draw', mapCoord, /Method, /AddCmd ENDIF IF Arg_Present(outimage) THEN BEGIN IF N_Elements(thisImage) NE 0 THEN BEGIN IF N_Elements(r) NE 0 THEN palette=[[r],[g],[b]] outimage = Temporary(thisImage) ENDIF ELSE BEGIN IF N_Elements(geofile) EQ 0 THEN BEGIN geofile = image void = Query_Tiff(geofile, info) IF info.channels GE 3 THEN BEGIN outimage = Read_Tiff(geofile, SUB_RECT=sub_rect) ENDIF ELSE BEGIN outimage = Read_Tiff(geofile, r, g, b, SUB_RECT=sub_rect) ENDELSE IF N_Elements(r) NE 0 THEN palette=[[r],[g],[b]] dims = Image_Dimensions(outimage, YINDEX=yindex) outimage = Reverse(Temporary(outimage), yindex+1) ENDIF ELSE outimage = image ENDELSE ENDIF ; Get return values. mapCoord -> GetProperty, $ BOUNDARY=boundary, $ ELLIPSOID=ellipsoid, $ LATLONBOX=latlonbox, $ MAP_PROJECTION=map_projection RETURN, mapCoord END