pro modis_map

filename='MYD021KM.A2015305.0030.006.2015305162839.hdf'

;  Open HDF file in SD mode
file_id=hdf_sd_start(filename,/read)

;  Get SD variable
var_id=hdf_sd_nametoindex(file_id,'Latitude')
sds_id=hdf_sd_select(file_id,var_id)
hdf_sd_getdata,sds_id,lat

var_id=hdf_sd_nametoindex(file_id,'Longitude')
sds_id=hdf_sd_select(file_id,var_id)
hdf_sd_getdata,sds_id,lon

var_id=hdf_sd_nametoindex(file_id,'EV_250_Aggr1km_RefSB')
sds_id=hdf_sd_select(file_id,var_id)
var_id=HDF_SD_ATTRFIND(sds_id,"reflectance_scales")
HDF_SD_ATTRINFO,sds_id,var_id,DATA=scale_factor
var_id=HDF_SD_ATTRFIND(sds_id,"reflectance_offsets")
HDF_SD_ATTRINFO,sds_id,var_id,DATA=offset


hdf_sd_getdata,sds_id,band_250

;  Close SD variable
hdf_sd_endaccess,sds_id

;  Close hdf file
hdf_sd_end,file_id

;dataValue = scale_factor * (dataValue-offset)

;  Get size of the data arrays
s=size(band_250,/dimensions)
sx_250=s[0] & sy_250=s[1]

;  Resize and linearly interpolate lat/lon arrays to data array (sds) size
lat_250=congrid(lat,sx_250,sy_250,/interp)
lon_250=smart_interp(lon,sx_250,sy_250)

;********Explore data

mydevice=!D.Name
pxdim=600
pydim=600
set_plot,'Z'
device,set_resolution=[pxdim,pydim]

;  Position plots
xl=0.13 & xr=0.85
yb=0.10 & yt=0.90
sx=0.05
sy=0.12
numplots_x=1
numplots_y=1
position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
cbpos=pos
cbpos[*,0]=cbpos[*,2]+0.08
cbpos[*,2]=cbpos[*,0]+0.02

pnum=0

d=convert_coord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
isx=(d[0,1,0]-d[0,0,0])
isy=(d[1,1,0]-d[1,0,0])
px0=d[0,0,0]
py0=d[1,0,0]

i1=image(lat,/buffer,$
    rgb_table=39,$
    dimensions=[pxdim,pydim],$  ;the size of the overall window
    image_dimensions=[isx,isy],$
    image_location=[px0,py0],$
    position=reform(pos[pnum,*]))

ct1=contour(lat,/current,$
    position=reform(pos[pnum,*]),$
    xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,$
    axis_style=2,$
    title='Latitude')

c0=colorbar(target=i1,$
    orientation=1,$
    position=[pos[pnum,2]+0.11,pos[pnum,1]+0.01,pos[pnum,2]+0.13,pos[pnum,3]-0.01],$
    tickdir=1,border_on=1)

i1.save,'lat.png',height=pydim,resolution=600

pnum=0

d=convert_coord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
isx=(d[0,1,0]-d[0,0,0])
isy=(d[1,1,0]-d[1,0,0])
px0=d[0,0,0]
py0=d[1,0,0]

i1=image(lon,/buffer,$
  rgb_table=39,$
  dimensions=[pxdim,pydim],$  ;the size of the overall window
  image_dimensions=[isx,isy],$
  image_location=[px0,py0],$
  position=reform(pos[pnum,*]))

ct1=contour(lon,/current,$
  position=reform(pos[pnum,*]),$
  xtickdir=1,ytickdir=1,$
  xstyle=1,ystyle=1,$
  axis_style=2,$
  title='Longitude')

c0=colorbar(target=i1,$
  orientation=1,$
  position=[pos[pnum,2]+0.11,pos[pnum,1]+0.01,pos[pnum,2]+0.13,pos[pnum,3]-0.01],$
  tickdir=1,border_on=1)

i1.save,'lon.png',height=pydim,resolution=600

pnum=0

d=convert_coord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
isx=(d[0,1,0]-d[0,0,0])
isy=(d[1,1,0]-d[1,0,0])
px0=d[0,0,0]
py0=d[1,0,0]

i1=image(band_250[*,*,0]*scale_factor[0],/buffer,$
  rgb_table=39,$
  dimensions=[pxdim,pydim],$  ;the size of the overall window
  image_dimensions=[isx,isy],$
  image_location=[px0,py0],$
  position=reform(pos[pnum,*]))

ct1=contour(band_250[*,*,0]*scale_factor[0],/current,$
  position=reform(pos[pnum,*]),$
  xtickdir=1,ytickdir=1,$
  xstyle=1,ystyle=1,$
  axis_style=2,$
  title='band 1 reflectance (620-670 nm) red',/nodata)

c0=colorbar(target=i1,$
  orientation=1,$
  position=[pos[pnum,2]+0.11,pos[pnum,1]+0.01,pos[pnum,2]+0.13,pos[pnum,3]-0.01],$
  tickdir=1,border_on=1)

i1.save,'band1.png',height=pydim,resolution=600

pnum=0

d=convert_coord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
isx=(d[0,1,0]-d[0,0,0])
isy=(d[1,1,0]-d[1,0,0])
px0=d[0,0,0]
py0=d[1,0,0]

i1=image(band_250[*,*,1]*scale_factor[1],/buffer,$
  rgb_table=39,$
  dimensions=[pxdim,pydim],$  ;the size of the overall window
  image_dimensions=[isx,isy],$
  image_location=[px0,py0],$
  position=reform(pos[pnum,*]))

ct1=contour(band_250[*,*,1]*scale_factor[1],/current,$
  position=reform(pos[pnum,*]),$
  xtickdir=1,ytickdir=1,$
  xstyle=1,ystyle=1,$
  axis_style=2,$
  title='band 2 reflectance (841-876 nm)',/nodata)

c0=colorbar(target=i1,$
  orientation=1,$
  position=[pos[pnum,2]+0.11,pos[pnum,1]+0.01,pos[pnum,2]+0.13,pos[pnum,3]-0.01],$
  tickdir=1,border_on=1)

i1.save,'band2.png',height=pydim,resolution=600

;*******************MAP*******************
set_plot,mydevice
band1=band_250[*,*,0]*scale_factor[0]

minlat=min(lat)
maxlat=max(lat)
minlon=min(lon)
maxlon=max(lon)

;  Calculate regular grid array
delta=0.05
ynum=(maxlat-minlat)/delta
grid_lat=findgen(ynum)*delta+minlat
xnum=(maxlon-minlon)/delta
grid_lon=findgen(xnum)*delta+minlon

;  Triangulate the data
triangulate,lon_250,lat_250,triangles,b

print,'Start gridding'
reg_grid=griddata(lon_250,lat_250,band1,/degree,$
  /grid,xout=grid_lon,yout=grid_lat,$
  method='NaturalNeighbor',missing=0,triangles=triangles)
  
pnum=0

i0=image(reg_grid,grid_lon,grid_lat,/buffer,$
  rgb_table=39,$
  dimensions=[pxdim,pydim],$
  grid_units='degrees',$
  limit=[minlat,minlon,maxlat,maxlon],$
  image_location=[minlon,minlat],$
  image_dimensions=[(maxlon-minlon),(maxlat-minlat)],$
  title='band 1 on a map',$
  irregular=0,$
  max_value=1.2,$
  min_value=0,$
  map_projection='mercator',$
  position=reform(pos[pnum,*]))
  
  i1=mapcontinents(color='white',/hires)
  i0.mapgrid.box_axes=1
  i0.mapgrid.label_position=0
  i0.mapgrid.box_color='black'
  i0.mapgrid.color='white'
  i0.mapgrid.label_color='black'
  i0.mapgrid.linestyle='dot'  

  c0=colorbar(target=i0,$
    orientation=1,$
    position=[pos[pnum,2]+0.11,pos[pnum,1]+0.01,pos[pnum,2]+0.13,pos[pnum,3]-0.01],$
    tickdir=1,border_on=1)

  
i0.save,'map_band1.png',height=pydim,resolution=600

stop


end