pro read_global_snow_ice
set_plot, 'z' 
!p.font=-1
;
;  sample code showing how to read global snow/ice maps and draw a .png picture
;
;  maps are in latlon projection at 1/25 degree grid step
;  4500 lines and 9000 pixels per line
;  upper left corener of the map is at 90N 180W
;  map orientation:     west -- 0 -- east
;  
;    
;
dimx=9000
dimy=4500

latmin=-90
latmax=90
lonmin=-180
lonmax=180

openr,1, 'gl_4km_snow_ice_map_2010340'
b=bytarr(dimx,dimy)
readu,1,b
close,1

;     Meaning of byte values:
;              0: WATER
;              1: LAND
;              2: SNOW
;              3: ICE
;            200: UNDETERMINED, LAND
;            210: UNDETERMINED, WATER
;
;
; generate a map (snow, ice, land, water and undetermined)
;
  nx1=9000
  ny1=4500

  device, font='Times Bold', /tt_font, set_character_size=[1,1],             $
       set_resolution=[nx1,ny1]
  tv, b, order=1

  result=tvrd(0,0)

  ir=bytarr(256)
  ig=bytarr(256)
  ib=bytarr(256)


  ir[0]=0
  ig[0]=0      ; water
  ib[0]=250

  ir[1]=143
  ig[1]=143         ; land
  ib[1]=64

  ir[2]=250
  ig[2]=250         ; snow 
  ib[2]=250

  ir[3]=250
  ig[3]=250         ; ICE 
  ib[3]=0

  ir[210]=100
  ig[210]=100        ; UNDETERMINED OVER WATER
  ib[210]=100

  ir[200]=150
  ig[200]=150       ; undetermined over land
  ib[200]=150


  fnm='test_map.png'
  write_png, fnm, result, ir, ig, ib
set_plot, 'x'
return
end