; docformat = 'rst' ; ; NAME: ; cgNCDFMap ; ; PURPOSE: ; The purpose of this function is to translate map projection information ; found in a netCDF file 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 map projection information ; found in a netCDF file into a map coordinate object (cgMap) that can be ; used to georeference images with a map data coordinate system. The Coyote ; Library is required. ; ; :Categories: ; Graphics, Map Projections ; ; :Examples: ; To create a map coordinate from a netCDF file from NSIDC:: ; file = 'C:\IDL\data\netCDF\jun_2004_amsr_e_bt.nc' ; mapCoord = cgNCDFMap(file, /Continents, /Grid) ; ; :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/ ; ; :File_comments: ; It is not possible to have a one-to-one mapping between every netCDF 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 netCDF 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 netCDF files as possible and I will take pains to do so if ; I know they are needed. ; :History: ; Converted from old ncdf_Coord program in Catalyst Library to run with Coyote Graphics routines. ; 9 November 2011. David W. Fanning. ; ; :Copyright: ; Copyright (c) 2011, Fanning Software Consulting, Inc. ;---------------------------------------------------------------------------------------- ; ;+-------------------------------------------------------------------------- ; :Description: ; Finds the correct projected map unit from the file. Assumes meters. ; ; :Params: ; varWithMap: in, required, type=string ; The name of the variable in the netCDF file that contains map ; projection information. ; fileObj: in, required, type=object ; The netCDF file object (ncdf_file) which can parse the netCDF file. ; ; :Keywords: ; silent: in, optional, type=boolean, default=0 ; If set, do not report errors. ;--------------------------------------------------------------------------- FUNCTION cgNCDFMap_GetMapUnit, varWithMap, fileObj, SILENT=silent ; 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() RETURN, mapUnit ENDIF ; Assume meters. mapUnit = 1. ; What kind of units are being used by the variable with the map projection? IF fileObj -> HasVarAttr(varWithMap, 'units') THEN BEGIN units = fileObj -> GetVarAttrValue(varWithMap, 'units') CASE StrLowCase(units) OF 'm': mapUnit = 1. 'metre': mapUnit = 1. 'meter': mapUnit = 1. 'meters': mapUnit = 1. 'km': mapUnit = 1000. 'kilometer': mapUnit = 1000. 'kilometers': mapUnit = 1000. ELSE: Message, 'Unrecognized unit of measurement: ' + units ENDCASE ENDIF RETURN, mapUnit END ;+-------------------------------------------------------------------------- ; :Description: ; Finds the boundary (XRANGE and YRANGE) of the map coordinate space ; and assigns these to the map object. ; ; :Params: ; varWithMap: in, required, type=string ; The name of the variable in the netCDF file that contains map ; projection information. ; thisMapVar: in, required, type=string ; The name of the map projection variable identified by the ; grid_mapping attribute of the varWithMap variable. ; fileObj: in, required, type=object ; The netCDF file object (ncdf_file) which can parse the netCDF file. ; mapCoord: in, required, type=object ; The map coordinate object we are creating. The XRange and YRange ; of the discovered boundaries are assigned to this object. ; ; :Keywords: ; silent: in, optional, type=boolean, default=0 ; If set, do not report errors. ; use_latlon: in, optional, type=boolean, default=0 ; If this keyword is set, the boundary ranges will be forced to be determined ; from the latitude and longitude arrays in the file. ; xrange: out, optional, type=double ; The discovered X projected meter range of the map projection. ; yrange: out, optional, type=double ; The discovered Y projected meter range of the map projection. ;--------------------------------------------------------------------------- FUNCTION cgNCDFMap_FindBoundary, varWithMap, thisMapVar, fileObj, mapCoord, $ SILENT=silent, XRANGE=xrange, YRANGE=yrange, USE_LATLON=use_latlon ; 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, success ENDIF ; Assume success success = 1 IF ~Keyword_Set(use_latlon) THEN BEGIN ; Let's look for helpful values in the projection variable. IF fileObj -> HasVarAttr(thisMapVar, 'grid_boundary_top_projected_y') THEN BEGIN y1 = fileObj -> GetVarAttrValue(thisMapVar, 'grid_boundary_top_projected_y') mapUnit = cgNCDFMap_GetMapUnit(thisMapVar, fileObj, SILENT=silent) y1 = y1 / mapUnit ENDIF IF fileObj -> HasVarAttr(thisMapVar, 'grid_boundary_bottom_projected_y') THEN BEGIN y0 = fileObj -> GetVarAttrValue(thisMapVar, 'grid_boundary_bottom_projected_y') mapUnit = cgNCDFMap_GetMapUnit(thisMapVar, fileObj, SILENT=silent) y0 = y0 / mapUnit ENDIF IF fileObj -> HasVarAttr(thisMapVar, 'grid_boundary_left_projected_x') THEN BEGIN x0 = fileObj -> GetVarAttrValue(thisMapVar, 'grid_boundary_left_projected_x') mapUnit = cgNCDFMap_GetMapUnit(thisMapVar, fileObj, SILENT=silent) x0 = x0 / mapUnit ENDIF IF fileObj -> HasVarAttr(thisMapVar, 'grid_boundary_right_projected_x') THEN BEGIN x1 = fileObj -> GetVarAttrValue(thisMapVar, 'grid_boundary_right_projected_x') mapUnit = cgNCDFMap_GetMapUnit(thisMapVar, fileObj, SILENT=silent) x1 = x1 / mapUnit ENDIF IF (Keyword_Set(x0) + Keyword_Set(x1) + Keyword_Set(y0) + Keyword_Set(y1)) EQ 4 THEN BEGIN xrange = [x0,x1] yrange = [y0,y1] mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange RETURN, success ENDIF ELSE Undefine, x0, x1, y0, y1 ENDIF ; Let's see if you can find a "latitude" and "longitude" variables. IF (fileObj -> HasVar('latitude', OBJECT=latVar)) AND $ (fileObj -> HasVar('longitude', OBJECT=lonVar)) THEN BEGIN lats = latVar -> GetValue() lons = lonVar -> GetValue() s = Size(lats, /DIMENSION) ; Convert to XY coords xy = Map_Proj_Forward(-180.0 > lons < 180.0, -90.0 > lats < 90.0, $ MAP_STRUCTURE=mapCoord->GetMapStructure()) x = Reform(xy[0,*], s[0], s[1]) y = Reform(xy[1,*], s[0], s[1]) ygrid = Reverse(Reform(y[0,*])) xgrid = Reform(x[*,0]) half_x_pixel = Abs(Median(xgrid - Shift(xgrid,1))) / 2.0 half_y_pixel = Abs(Median(ygrid - Shift(ygrid,1))) / 2.0 xrange = [xgrid[0]-half_x_pixel, xgrid[s[0]-1]+half_x_pixel] yrange = [ygrid[0]-half_y_pixel, ygrid[s[1]-1]+half_y_pixel] mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange RETURN, success ENDIF ; Can we stop now? IF (N_Elements(xrange) NE 0) AND (N_Elements(yrange) NE 0) THEN BEGIN mapCoord -> SetProperty, XRANGE=xrange, YRANGE=yrange RETURN, success ENDIF ; Look in the data variable for the attribute "coordinates" IF fileObj -> HasVarAttr(varWithMap, 'coordinates', OBJECT=coordAttr) THEN BEGIN coordVars = coordAttr -> GetValue() ; This should be the name of two variables. parts = StrSplit(coordVars, /Extract) IF N_Elements(parts) NE 2 THEN Message, 'Coordinates attribute did not yield two variable names.' FOR j=0,1 DO BEGIN thisCoordVar = parts[j] IF fileObj -> HasVarAttr(thisCoordVar, 'standard_name') THEN BEGIN thisCoordName = fileObj -> GetVarAttrValue(thisCoordVar, 'standard_name') ENDIF ELSE thisCoordName = "No Standard Name" IF fileObj -> HasVarAttr(thisCoordVar, 'units') THEN BEGIN mapUnit = cgNCDFMap_GetMapUnit(thisCoordVar, fileObj, SILENT=silent) ENDIF ; Hopefully, this is a projection coordinate. CASE thisCoordName OF 'projection_x_coordinate': BEGIN xvec = fileObj -> GetVarData(thisCoordVar) nvec = N_Elements(xvec) spacing = Median(xvec - Shift(xvec, 1)) halfpixel = spacing/2.0 xrange = [xvec[0]-halfpixel, xvec[nvec-1]+halfpixel] / mapUnit END 'projection_y_coordinate': BEGIN yvec = fileObj -> GetVarData(thisCoordVar) nvec = N_Elements(yvec) spacing = Median(yvec - Shift(yvec, 1)) halfpixel = spacing/2.0 yrange = [yvec[0]-halfpixel, yvec[nvec-1]+halfpixel] / mapUnit END ELSE: Message, 'Do not recognize the coordinate name ' + thisCoordName + '.' ENDCASE ENDFOR ENDIF ELSE Return, 0 ; Didn't find anything. END ;+-------------------------------------------------------------------------- ; :Description: ; Finds the length of the semi-major and semi-minor axes of the ; ellipsoid used in the netCDF file. ; ; :Params: ; varWithMap: in, required, type=string ; The name of the variable in the netCDF file that contains map ; projection information. ; thisMapVar: in, required, type=string ; The name of the map projection variable identified by the ; grid_mapping attribute of the varWithMap variable. ; fileObj: in, required, type=object ; The netCDF file object (ncdf_file) which can parse the netCDF file. ; ; :Keywords: ; silent: in, optional, type=boolean, default=0 ; If set, do not report errors. ; semimajor_axis: out, optional, type=double ; The length of the semi-major axis. ; semiminor_axis: out, optional, type=double ; The length of the semi-minor axis. ;--------------------------------------------------------------------------- FUNCTION cgNCDFMap_EllipseAxes, varWithMap, thisMapVar, fileObj, SILENT=silent, $ SEMIMAJOR_AXIS=semimajor_axis, SEMIMINOR_AXIS=semiminor_axis ; varWithMap -- The name of the variable containing the grid_mapping attribute. ; thisMapVar -- The map projection variable identified by the grid_mapping attribute. ; fileObj -- The NCDF_FILE object for the file. ; ; SILENT -- If set, errors are handled silently. ; SEMIMAJOR_AXIS -- Output, the semi-major axis length in meters. ; SEMIMINOR_AXIS -- Output, the semi-major axis length in meters. ; 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, success ENDIF ; Assume success success = 1 ; Convert whatever map units are associated with this variable ; into a mapUnit that I can use to divide lengths by to get meters. mapUnit = cgNCDFMap_GetMapUnit(varWithMap, fileObj, SILENT=silent) ; Look for an "earth_radius" attribute. IF fileObj -> HasVarAttr(thisMapVar, 'earth_radius') THEN BEGIN radius = fileObj -> GetVarAttrValue(thisMapVar, 'earth_radius') semimajor_axis = radius/mapUnit semiminor_axis = radius/mapUnit ENDIF ; Are you done? IF (Keyword_Set(semimajor_axis) AND Keyword_Set(semiminor_axis)) EQ 2 THEN RETURN, 1 ; Try looking for CF 1.4 compliant attributes "semi_major_axis" and "semi_minor_axis". IF fileObj -> HasVarAttr(thisMapVar, 'semi_major_axis') THEN BEGIN semimajor_axis = fileObj -> GetVarAttrValue(thisMapVar, 'semi_major_axis') semimajor_axis = semimajor_axis / mapUnit ENDIF IF fileObj -> HasVarAttr(thisMapVar, 'semi_minor_axis') THEN BEGIN semiminor_axis = fileObj -> GetVarAttrValue(thisMapVar, 'semi_minor_axis') semiminor_axis = semiminor_axis / mapUnit ENDIF ELSE BEGIN IF N_Elements(semimajor_axis) NE 0 THEN BEGIN IF fileObj -> HasVarAttr(thisMapVar, 'inverse_flattening') THEN BEGIN f_inv = fileObj -> GetVarAttrValue(thisMapVar, 'inverse_flattening') f = 1.0D/f_inv semiminor_axis = semimajor_axis - ( f * semimajor_axis ) ENDIF ; Perhaps this is a sphere. IF N_Elements(semiminor_axis) EQ 0 THEN semiminor_axis = semimajor_axis ENDIF ENDELSE ; Are you done? IF (Keyword_Set(semimajor_axis) AND Keyword_Set(semiminor_axis)) EQ 2 THEN RETURN, success ; Try looking for CF 1.4 non-compliant attributes. IF fileObj -> HasVarAttr(thisMapVar, 'semimajor_radius') THEN BEGIN semimajor_axis = fileObj -> GetVarAttrValue(thisMapVar, 'semimajor_radius') ENDIF IF fileObj -> HasVarAttr(thisMapVar, 'semiminor_radius') THEN BEGIN semiminor_axis = fileObj -> GetVarAttrValue(thisMapVar, 'semiminor_radius') ENDIF ELSE BEGIN IF N_Elements(semimajor_axis) NE 0 THEN BEGIN IF fileObj -> HasVarAttr(thisMapVar, 'inverse_flattening') THEN BEGIN f_inv = fileObj -> GetVarAttrValue(thisMapVar, 'inverse_flattening') f = 1.0D/f_inv semiminor_axis = semimajor_axis - ( f * semimajor_axis ) ENDIF ; Perhaps this is a sphere. IF N_Elements(semiminor_axis) EQ 0 THEN semiminor_axis = semimajor_axis ENDIF ENDELSE ; Are you done? IF (Keyword_Set(semimajor_axis) AND Keyword_Set(semiminor_axis)) EQ 2 THEN RETURN, success ; Maybe there is an ESRI Projection Engine String we can use. IF fileObj -> HasVarAttr(varWithMap, 'esri_pe_string') THEN BEGIN esri_str = fileObj -> GetVarAttrValue(varWithMap, 'esri_pe_string') pos = StrPos(esri_str, 'SPHEROID[') IF pos NE -1 THEN BEGIN substring = StrMid(esri_str, pos) closebracket = StrPos(substring, ']') spheroid = StrMid(substring, 9, closebracket-9) parts = StrSplit(spheroid, ',', /EXTRACT) semimajor_axis = Double(parts[1]) IF mapUnit EQ 1000 THEN semimajor_axis = semimajor_axis / mapUnit f_inv = Double(parts[2]) f = 1.0D/f_inv semiminor_axis = semimajor_axis - ( f * semimajor_axis ) ENDIF ENDIF ; Did we find anything? IF (Keyword_Set(semimajor_axis)EQ 0) OR $ (Keyword_Set(semiminor_axis) EQ 0) THEN BEGIN ; Hells bells! Let's just issue a warning and assume WGS-84. What else ; can we do!? Message, /INFORMATIONAL, 'Cannot find any sensible datum radius in file. Assuming WGS-84 values.' semimajor_axis = 6378137.0D semiminor_axis = 6356752.31414D RETURN, success ENDIF ELSE RETURN, success END ;+--------------------------------------------------------------------------------------- ; :Description: ; The purpose of this function is to translate map projection information ; found in a netCDF file into a map coordinate object (cgMap) that can be ; used to georeference images with a map data coordinate system. The Coyote ; Library is required. ; ; :Params: ; ncdf_filename: in, optional, type=string ; The name of a netCDF file containing map projection information from ; which a cgMap object can be created. ; ; :Keywords: ; 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. ; 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. ; 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. ; 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. ; 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. ; 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!) ; 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. ; use_latlon: in, optional, type=boolean, default=0 ; If this keyword is set, the boundary ranges will be forced to be determined ; from the latitude and longitude arrays in the file. ; xrange: out, optional, type=double ; The discovered X projected meter range of the map projection. ; yrange: out, optional, type=double ; The discovered Y projected meter range of the map projection. ;------------------------------------------------------------------------------------ FUNCTION cgNCDFMap, ncdf_filename, $ GCOLOR=gcolor, $ GRID=grid, $ CONTINENTS=continents, $ CCOLOR=ccolor, $ MCOLOR=mcolor, $ ONIMAGE=onimage, $ SILENT=silent, $ SUCCESS=success, $ TITLE=title, $ USE_LATLON=use_latlon, $ XRANGE=xrange, $ YRANGE=yrange ; Compiler options. Compile_Opt DEFINT32 Compile_Opt STRICTARR Compile_Opt STRICTARRSUBS Compile_Opt LOGICAL_PREDICATE Catch, theError IF theError NE 0 THEN BEGIN Catch, theError IF ~Keyword_Set(silent) THEN void = cgErrorMsg() success = 0 IF Obj_Valid(fileObj) THEN Obj_Destroy, fileObj RETURN, Obj_New() ENDIF ; Need a filename? IF N_Elements(ncdf_filename) EQ 0 THEN BEGIN ncdf_filename = Dialog_Pickfile(TITLE='Select netCDF File...', FILTER=['*.nc','*.ncdf']) IF ncdf_filename EQ "" THEN RETURN, Obj_New() ENDIF ; Set default values of keywords. SetDefaultValue, gcolor, 'gray' SetDefaultValue, ccolor, 'charcoal' SetDefaultValue, mcolor, 'black' ; Assume success success = 1 ; Create a NCDF_FILE object for reading the data file. fileObj = Obj_New('NCDF_FILE', ncdf_filename, /NOCLUTTER) IF Obj_Valid(fileObj) EQ 0 THEN Message, 'Cannot create NCDF_FILE object for reading data file.' ; Get a list of the variables in the file. Check each list. The first ; one that has a "grid_mapping" attribute will be the map projection we ; use here. varNames = fileObj -> GetVarNames(COUNT=varCount) FOR j=0,varCount-1 DO BEGIN IF ~fileObj->HasVarAttr(varNames[j], 'grid_mapping') THEN Continue varWithMap = varNames[j] thisMapVar = fileObj -> GetVarAttrValue(varWithMap, 'grid_mapping') thisMapProj = fileObj -> GetVarAttrValue(thisMapVar, 'grid_mapping_name') IF thisMapProj EQ "" THEN Message, 'Cannot find map projection name in file.' BREAK ENDFOR IF N_Elements(varWithMap) EQ 0 THEN $ Message, 'This netCDF file has no variables with map projection information.' ; Define map projection variables you will find in netCDF files. cf_std_names = [{name:'albers_conical_equal_area', possible:1, gctp:103}, $ {name:'azimuthal_equidistant', possible:1, gctp:112}, $ {name:'lambert_azimuthal_equal_area', possible:1, gctp:111}, $ {name:'lambert_conformal_conic', possible:1, gctp:104}, $ {name:'lambert_cylindrical_equal_area', possible:0, gctp:-1}, $ {name:'latitude-longitude', possible:0, gctp:-1}, $ {name:'mercator', possible:1, gctp:105}, $ {name:'orthographic', possible:1, gctp:114}, $ {name:'polar_stereographic', possible:1, gctp:106}, $ {name:'rotated_latitude_longitude', possible:0, gctp:-1}, $ {name:'stereographic', possible:1, gctp:110}, $ {name:'transverse _mercator', possible:1, gctp:109}, $ {name:'vertical_perspective', possible:0, gctp:-1}] ; Define the expected parameters for those netCDF projection types. cf_std_name_parameters = [ Ptr_New(['standard parallel', $ ; Albers Equal Area 'longitude_of_central_meridian', $ 'latitude_of_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['longitude_of_projection_origin', $ ; Azimuthal Equidistant 'latitude_of_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['longitude_of_projection_origin', $ ; Lambert Azimutal Equal Area 'latitude_of_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['standard parallel', $ ; Lambert Conformal 'longitude_of_projection_origin', $ 'latitude_of_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['longitude_of_central_meridian', $ ; Lambert Cylindrical Equal Area 'standard parallel', $ 'scale_factor_at_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(), $ ; Latitude-Longitude Ptr_New(['longitude_of_projection_origin', $ ; Mercator 'standard parallel', $ 'latitude_of_true_scale', $ 'scale_factor_at_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['longitude_of_projection_origin', $ ; Orthographic 'latitude_of_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['straight_vertical_longitude_from_pole', $ ; Polar Stereographic 'longitude_of_projection_origin', $ 'latitude_of_projection_origin', $ 'standard parallel', $ 'latitude_of_true_scale', $ 'scale_factor_at_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(), $ ; Rotated Pole Ptr_New(['longitude_of_projection_origin', $ ; Stereographic 'latitude_of_projection_origin', $ 'scale_factor_at_projection_origin', $ 'false_easting', 'false_northing']), $ Ptr_New(['scale_factor_at_central_meridian', $ ; Transverse Mercator 'longitude_of_central_meridian', $ 'latitude_of_central_meridian', $ 'false_easting', 'false_northing']), $ Ptr_New() ] ; Vertical Perspective ; Can you locate the map projection in the list of projections you know about? index = Where(cf_std_names.name EQ thisMapProj, count) IF count EQ 0 THEN Message, 'Cannot process map projection: ' + thisMapProj ; Is is possible to do this map projection in IDL? IF cf_std_names[index].possible EQ 0 THEN $ Message, 'Cannot process map projection "' + thisMapProj + '" in IDL.' ; Can you find values for the ellipsoid? success = cgNCDFMap_EllipseAxes(varWithMap, thisMapVar, fileObj, SILENT=silent, $ SEMIMAJOR_AXIS=semimajor_axis, SEMIMINOR_AXIS=semiminor_axis) IF ~success THEN Message, 'Cannot locate ellipsoid axes in file.' ; Get the attributes for this particular map projection. attrs = *(cf_std_name_parameters[index])[0] nAttrs = N_Elements(attrs) FOR j=0,nAttrs-1 DO BEGIN thisAttr = attrs[j] IF fileObj -> HasVarAttr(thisMapVar, thisAttr, OBJECT=attrObj) THEN BEGIN attrValue = attrObj -> GetValue() ENDIF ELSE Continue CASE thisAttr OF 'longitude_of_central_meridian' : center_longitude = attrValue 'latitude_of_central_meridian' : center_latitude = attrValue 'false_easting' : false_easting = attrValue 'false_northing' : false_northing = attrValue 'longitude_of_central_meridian' : center_longitude = attrValue 'latitude_of_projection_origin' : center_latitude = attrValue 'longitude_of_projection_origin' : center_longitude = attrValue 'latitude_of_true_scale': BEGIN CASE cf_std_names[index].name OF 'polar_stereographic': center_latitude = attrValue ELSE: true_scale_latitude = attrValue ENDCASE END 'straight_vertical_longitude_from_pole': center_longitude = attrValue 'standard_parallel' : BEGIN IF N_Elements(attrValue) EQ 2 THEN BEGIN standard_par1 = attrValue[0] standard_par2 = attrValue[1] ENDIF ELSE standard_par1 = attrValue END 'scale_factor_at_projection_origin': BEGIN CASE cf_std_names[index].name OF 'transverse _mercator': mercator_scale=attrValue 'mercator': mercator_scale=attrValue ELSE: Message, /Informational, 'Not sure what to do with the attribute ' + $ '"scale_factor_at_projection_origin".' ENDCASE END 'scale_factor_at_central_meridian': BEGIN CASE cf_std_names[index].name OF 'transverse _mercator': mercator_scale=attrValue 'mercator': mercator_scale=attrValue ELSE: Message, /Informational, 'Not sure what to do with the attribute ' + $ '"scale_factor_at_central_meridian".' ENDCASE END ENDCASE ENDFOR mapCoord = Obj_New('cgMap', cf_std_names[index].gctp, $ CENTER_LATITUDE=center_latitude, $ CENTER_LONGITUDE=center_longitude, $ COLOR=mcolor, $ SEMIMAJOR_AXIS=semimajor_axis, $ SEMIMINOR_AXIS=semiminor_axis, $ NAME=name, $ MERCATOR_SCALE=mercator_scale, $ ONIMAGE=Keyword_Set(onimage), $ STANDARD_PAR1=standard_par1, $ STANDARD_PAR2=standard_par2, $ TITLE=title, $ TRUE_SCALE_LATITUDE=true_scale_latitude, $ FALSE_EASTING=false_easting, $ FALSE_NORTHING=false_northing) ; Now, can we find the extend of the XY Cartesian boundaries to set the range ; of the coordinate object. success = cgNCDFMap_FindBoundary(varWithMap, thisMapVar, fileObj, mapCoord, $ SILENT=silent, USE_LATLON=use_latlon, XRANGE=xrange, YRANGE=yrange) IF ~success THEN Message, 'Cannot find XY boundary values in file.' ; 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 ; Destroy the file object you created. Obj_Destroy, fileObj RETURN, mapCoord END