Internal: Example how to display full-disk Goes-r image using IDL

Description:

This article provides a demonstration how to use the IMAGE function to display a full-disk Goes-r image. 

Solution:

The procedure to display a full-disk Goes-r image using the IMAGE function is shown below:

1) Read the "Rad" data out of the file.
2) Read the "X" and "Y" values from the file and converted them to Radians
3) Use the set of equations from a Goes-R spec that I found online to convert the X and Y data to Latitude and Longitude values.
4) Convert the Lat/Lon coordinates calculated in step 3 to UV coordinates (Cartesian).
5) Use the "Rad" data and the  UV coordinates to display the image. 
 
An example set of code that shows how this can be achieved is shown below:

function dj_goes_attempt_2_sx, rs, x, y
  compile_opt idl2
  
  sx =  rs*cos(x)*cos(y)
  return, sx
  
end


function dj_goes_attempt_2_sy, rs, x, y
  compile_opt idl2

  sy = rs*sin(x)
  return, sy
  
end

function dj_goes_attempt_2_sz, rs, x, y
  compile_opt idl2

  sz = rs*cos(x)*sin(y)
  return, sz
  
end

function dj_goes_attempt_2_lat, req, rpol, H, sx, sy, sz
    compile_opt idl2

    num = (req^2)*sz
    den = (rpol^2)*(sqrt((H-sx)^2 + sy^2))
    lat = ATAN(num/den)
    return, lat    
    
end

function dj_goes_attempt_2_lon, H, sx, sy
  compile_opt idl2

   lon = atan(sy/(H-sx))  
   return, lon
  
end



function dj_goes_attempt_2_rs, b, a, c
  compile_opt idl2

   rs = (-b-sqrt(b^2 - 4*a*c))/(2*a)
   return, rs
   
end


function dj_goes_attempt_2_c, H, req
  compile_opt idl2

  c = H^2 - req^2
  return, c
  
end


function dj_goes_attempt_2_b, H, x, y
  compile_opt idl2

   b = -2*H*cos(x)*cos(y)
   return, b

end


function dj_goes_attempt_2_a, x, y, req, rpol
  compile_opt idl2

   a = sin(x)^2 + (cos(x)^2) * (cos(y)^2 + ((req^2)/(rpol^2))*sin(y)^2)
  return, a
  
end

; Live your life in fea
pro dj_goes_attempt_2
  compile_opt idl2

  ;Open the Goes-R file
  fw = file_which('OR_ABI-L1b-RadF-M4C03_G16_s20153352323571_e20153352328390_c20153352328450.nc')
  
  
  ;Read in the "Rad" data 
  ;NOTE: my computer not very hardcore so I 
  ;           had to sample data set to get the processing 
  ;           to work within a decent amount of time 
  NCDF_GET, fw, "Rad", result
  rad_hash = result['Rad']
  rad_values = rad_hash['value'] 
  rad_values = congrid(rad_values, (10848/5), (10848/5)) 
  
  ;Get the X Radians data
  ncdf_get,fw,"x",x_res
  x_hash = x_res['x']
  x_array = x_hash['value']
  x_rads = x_array*(2.80000e-005) - 0.151858
  x_rads = x_rads[0:10847:5]
  x_rads = transpose(congrid(transpose(x_rads), 10848/5, 10848/5))
  x_rads = congrid(x_rads, 10848/5, 10848/5)

 ;Get the Y Radians data
  ncdf_get,fw,'y',y_res
  y_hash = y_res['y']
  y_array = y_hash['value']
  y_rads = y_array*(-2.80000e-005) + 0.151858
  y_rads = y_rads[0:10847:5]
  y_rads = congrid(TRANSPOSE(y_rads), 10848/5, 10848/5)
  y_rads = congrid(y_rads, 10848/5, 10848/5)

  ;Define a few variables that can be found
  ;in the file
  req = 6378137.0
  rpol = 6356752.3
  h = 42164160. ;35786023.
    
  ;Use the equations from the specification document to 
  ;convert the X/Y Radian data to Lon/Lat data
  a = dj_goes_attempt_2_a(x_rads, y_rads, req, rpol)
  b =  dj_goes_attempt_2_b(h, x_rads, y_rads)
  c = dj_goes_attempt_2_c( h, req)
  
  rs = dj_goes_attempt_2_rs(b, a, c)
  sx = dj_goes_attempt_2_sx(rs, x_rads, y_rads)
  sy = dj_goes_attempt_2_sy(rs, x_rads, y_rads)
  sz = dj_goes_attempt_2_sz(rs, x_rads, y_rads)

  lat = (dj_goes_attempt_2_lat(req, rpol, H, sx, sy, sz))*!RADEG 
  lon =  dj_goes_attempt_2_lon(H, sx, sy)*!RADEG-75
  
  ;Covert lat/lon data to UV coordinates (catresian). I think that IDL
  ;should be able to do this under the covers but for some reason it 
  ;seems to bonk and then deform the image
  me_map = map('near side perspective', CENTER_LONGITUDE=-75, /HIDE, /buffer)
  map_xy = me_map.mapforward(reform(lon), reform(lat)) 
  
  xx = map_xy[0,*]
  yy = map_xy[1,*]
  
  xx = reform(xx, 2169, 2169)
  yy = reform(yy, 2169, 2169)
  
  ;Display the Image using GRID_UNITS=1. The country lines appear to match pretty well. 
   i = image(reverse(rad_values,2), xx, yy, MAP_PROJECTION='near side perspective', max_value=1022, MIN_VALUE=0 , GRID_UNITS=1,CENTER_LONGITUDE=-75)
  mc = MAPCONTINENTS(/COUNTRIES, color='light blue', thick=3, /hires)
  

end