; Antarctic 1-km Digital Elevation Model (DEM) from 
; Combined ERS-1 Radar and ICESat Laser Satellite Altimetry
;
;  ftp://sidads.colorado.edu/pub/DATASETS/DEM/nsidc0422_antarctic_1km_dem/
;---
; ENVI:
;  map info = {Polar Stereographic South,
;              1.0000, 1.0000, -2800500.0000, 2800500.0000,
;              1.0000000000e+03, 1.0000000000e+03, WGS-84, units=Meters}
;  projection info = {31, 6378137.0, 6356752.3, -71.000000, 0.000000,
;              0.0, 0.0, WGS-84, Polar Stereographic South, units=Meters}
;
;---
;  http://nsidc.org/data/docs/daac/nsidc0422_antarctic_1km_dem/index.html
;---
; Time stamp: The DEM has a time stamp of January 2004.
;---
; Projection: World Geodetic System 1984 (WGS 84) is the reference ellipsoid and datum.
;
; The projection is polar stereographic with a standard parallel of 71S and 
; origin of the South Pole.
;---
; Grid Description: The grid is a 1 km grid of 5601 by 5601 grid boxes 
; centered symmetrically about the pole.
;---
; Elevation values in the DEM range from -999 to 4185.72 meters. 
; Values of -999 indicate grid cells with no measured or interpolated value, 
; and typically occur in areas beyond the geographic extent of Antarctic
;=============================================================================
;                            MAIN
;=============================================================================
 N    = 5601                                           ; see documentation

 diri = "./"
 fili = "krigged_dem_nsidc.bin"
 pthi = diri+fili

 setfileoption("bin","ReadByteOrder","LittleEndian")   ; see documentation
 
 dem   = fbindirread(pthi, 0, (/N,N/), "float")
 dem@_FillValue = -999.0                               ; see documentation

;---Add meta
 dem@long_name  = "Antarctic 1-km Digital Elevation Model (DEM)"
 dem@units  = "m"
 dem@info_1 = "Source: Combined ERS-1 Radar and ICESat Laser Satellite Altimetry"
 dem@info_2 = "Elevation in meters relative to WGS 84"

 printVarSummary(dem)
 printMinMax(dem,0)

;---Explore data
;opt = True
;opt@PrintStat = True
;statb = stat_dispersion(dem, opt )

;---GRAPHICS

 wks  = gsn_open_wks ("png", "DEM")                  ; send graphics to PNG file

 setvalues NhlGetWorkspaceObjectId() 
    "wsMaximumSize" : 200000000
 end setvalues 

 res                      = True                 ; plot mods desired
 res@cnFillOn             = True                 ; turn on color
 res@cnLinesOn            = False                ; no contour lines
 res@cnFillPalette        = "BlAqGrYeOrReVi200"  ; set color map
 res@cnFillMode           = "RasterFill"
 res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
 res@cnMinLevelValF       =    0.              ; set min contour level
 res@cnMaxLevelValF       = 4000.              ; set max contour level
 res@cnLevelSpacingF      =  250.              ; set contour spacing
 plot = gsn_csm_contour(wks,dem,res)