; docformat = 'rst' ;+ ; This is an example program to illustrate how to browse an HDF-EOS grid ; data file and display one of the data fields as a contour plot on a ; map projection grid from the file. ; ; :Params: ; ; file: in, required, type=string ; The input data file. For this example, the file 'AMSR_E_L3_RainGrid_B05_200707.hdf'. ;- PRO EOS_Grid_Display, file ; Need a file? IF N_Elements(file) EQ 0 THEN BEGIN file = cgPickfile(File='AMSR_E_L3_RainGrid_B05_200707.hdf') IF file EQ "" THEN RETURN ENDIF ; Determine how many grids are in the file, and what the grid names are. gridCount = EOS_GD_InqGrid(file, list) IF gridCount LE 0 THEN Message, 'This may not be a grid file. No grids found in file.' IF gridCount GT 1 THEN BEGIN grid = StrSplit(list, ',', /Extract) ENDIF ELSE grid = list Print, '' Print, 'The following grids are in the file: ' FOR j=0,N_Elements(grid)-1 DO Print, ' "' + grid[j] + '"' Print, '' ; Open the file for reading. fileID = EOS_GD_Open(file) ; Browse each grid in the file. FOR g=0,gridCount-1 DO BEGIN ; Make an attachment to the first grid. thisGrid = grid[g] gridID = EOS_GD_Attach(fileID, thisGrid) ; Determine the number and names of the attributes in this grid. attrCount = EOS_GD_InqAttrs(gridID, attrlist) IF attrCount LE 0 THEN Print, 'No attributes found in grid ' + thisGrid IF attrCount GT 1 THEN BEGIN attr = StrSplit(attrlist, ',', /Extract) ENDIF ELSE attr = attrlist Print, 'The following attributes are in the grid "' + thisGrid + '": ' FOR j=0,N_Elements(attr)-1 DO BEGIN Print, ' "' + attr[j] + '"' thisAttr = attr[j] ; This function is suppose to return a 0 or 1 to indicate success. ; But it ALWAYS returns a 0! success = EOS_GD_ReadAttr(gridID, thisAttr, value) IF success EQ 0 THEN Print, ' ', value[0] ENDFOR ; Get the projection information out of the grid. You MUST specify all these ; parameters, even if you don't intend to use them. void = EOS_GD_ProjInfo(gridID, projcode, zonecode, ellipsecode, projparams) projcodesIndex = [ 0, 1, 4, 6, 7, 9, 11, 20, 22, 24, 99] projcodesNames = ['Geographic', 'UTM', 'Lambert Conformal Conic', 'Polar Stereographic', $ 'Polyconic', 'Transverse Mercator', 'Lambert Equal Area', 'Hotine Oblique Mercator', $ 'Space Oblique Mercator', 'Goode Homolosine', 'Intergerized Sinusoidal'] index = Where(projcodesindex EQ projcode) thisProjection = projcodesNames[index] IF thisProjection EQ 'Intergerized Sinusoidal' THEN projcode = 31 IF thisProjection EQ 'Geographic' THEN projcode = 17 projcode = projcode + 100 ; To conform with Map_Proj_Init codes. Print, '' Print, ' Projection Information: ' Print, ' GCTP Projection: ', thisProjection Print, ' CGTP Proj Code: ', projcode Print, ' GCTP Zone: ', zonecode Print, ' GCTP Ellipsoid: ', ellipsecode ; Get the dimensions of the grid, and the locations of the upper-left and lower-right ; corners in latitude/longitude space. The values here are packed latitude/longitudes ; in DMS format. They must be divided by 10^6 to produce the correct value. void = EOS_GD_GridInfo(gridID, xsize, ysize, upleft, lowright) upleft = upleft / 10.^6 lowright = lowright / 10.^6 Print, '' Print, ' Dimensions of Grid: ', xsize, ysize Print, ' Upper-left Corner: ', upleft Print, ' Lower-right Corner: ', lowright ; Determine the number and names of the data fields in this grid. fieldCount = EOS_GD_InqFields(gridID, fieldlist, rank, type) IF fieldCount LE 0 THEN Message, 'No fields found in grid ' + thisGrid IF fieldCount GT 1 THEN BEGIN field = StrSplit(fieldlist, ',', /Extract) ENDIF ELSE field = fieldlist Print, '' Print, 'The following fields are in grid "' + thisGrid + '": ' FOR j=0,N_Elements(field)-1 DO Print, ' "' + field[j] + '"' ENDFOR ; Now that we know something about the file, we can visualize the fields in the file. ; Start with the field, "RrLandRain". Get the field data and be sure to ; reverse it, since we are working with map projections and IDL is backwards to ; the rest of the world. thisField = field[1] void = EOS_GD_ReadField(gridID, thisField, fieldValue) fieldValue = Reverse(fieldValue, 2) ; Try to find the missing value in the file. Clearly, this file does ;not include this information. (The success flag returns FAILED.) success = EOS_GD_GetFillValue(gridID, thisField, fillValue) Print, 'Fill value in file for field "' + thisField + '": ', fillValue Print, 'Success at locating Fill Value in file: ', (success EQ 0) indices = Where(fieldValue LT 0.0, cnt) IF cnt GT 0 THEN fieldValue[indices] = !Values.F_NaN ; Create vectors of latitudes and longitudes. lons = Scale_Vector(Findgen(xsize), upleft[0], lowright[0]) lats = Scale_Vector(Findgen(ysize), lowright[1], upleft[1]) ; Load the colors. nlevels = 12 cgLoadCT, 2, NColors=nlevels-1, Bottom=1, /Brewer, /Reverse TVLCT, cgColor('white', /Triple), nlevels ; Create a map projection space. position = [0.05, 0.075, 0.95, 0.78] cgDisplay, 800, 625 map = Obj_New('cgmap', projCode, Ellipsoid=ellipsoidCode, XRange=[-180, 180], $ YRange=[-80,80], Position=position, /LatLon_Range, $ Center_Longitude=0.0, Center_Latitude=0.0, Background='light gray') map -> Draw, /Erase ; Get the contour levels. Get the step, as you will need it below. levels = cgConLevels(fieldValue, NLEVELS=nlevels+1, MinValue=0.0, STEP=step) ; Draw the contour, taking care to pass the map object, so the latitude ; and longitude grids can be converted to projected meter space. cgContour, fieldValue, lons, lats, LEVELS=levels, $ C_Colors=Indgen(nlevels)+1, Position=position, $ XStyle=5, YStyle=5, /Cell_Fill, XTitle='Longitude', $ YTitle='Latitude', /Overplot, Map_Object=map ; Add a colorbar. cgColorbar, NColors=nlevels-1, Bottom=1, OOB_High=nlevels[0], $ Range=[Min(levels), Max(levels)-step], Divisions=nlevels[0],$ /Discrete, Title='Rain (mm)', Position=[0.88, 0.1, 0.92, 0.9] ; Add map annotations. cgMap_Continents, Color='tg1', Map=map cgMap_Grid, MAP=map, /Box cgText, 0.05, 0.025, /Normal, 'July 2007' ; Clean up Obj_Destroy, map void = EOS_GD_Detach(gridID) void = EOS_GD_Close(fileID) END