FUNCTION GRID_IMAGE, GRID, MISSING=MISSING, LONGITUDE=LONGITUDE

;- Convert a variable defined on the 25x25 km grid (dimension 814880)
;- to an image on a sinusoidal projection (dimensions 1600 x 800)

compile_opt idl2

;- Check arguments
if (n_elements(grid) eq 0) then message, 'Argument GRID is undefined'
dims = size(grid, /dimensions)
if (dims[0] ne 814880) then $
  message, 'Argument GRID must have dimensions [814880 elements]'
if (n_elements(missing) eq 0) then missing = -9999.0
if (n_elements(longitude) eq 0) then longitude = 0.0

;- Get first column from grid index file
openr, lun, 'get_grid_index.dat', /get_lun
header = ''
readf, lun, header, format='(a)'
data = fltarr(2, 800)
readf, lun, data
free_lun, lun
sum = reform(long(data[0, *]))

;- Create row map array
row = lonarr(814880)
loc1 = 0
for index = 0, 799 do begin
  loc2 = sum[index] - 1
  row[loc1 : loc2] = index
  loc1 = sum[index]
endfor

;- Create column map array
col = lonarr(814880)
loc1 = 0
for index = 0, 799 do begin
  loc2 = sum[index] - 1
  ncol = loc2 - loc1 + 1
  col[loc1 : loc2] = lindgen(ncol) + ((1600 - ncol) / 2)
  nshift = long(ncol * ((longitude + 180.0) / 360.0))
  col[loc1 : loc2] = shift(col[loc1 : loc2], nshift)
  loc1 = sum[index]
endfor

;- Transform the grid array to an image
image = replicate(missing, 1600, 800)
index = lindgen(814880)
image[col[index], row[index]] = grid
return, image

END