; Copyright (C) 2013 The HDF Group
;   All rights reserved.
;
;  This example code illustrates how to access and visualize 
; GES DISC OMI HDF-EOS5 Swath file with Google Earth using NCL. 
;
;  If you have any questions, suggestions, comments on this example, 
; please use the HDF-EOS Forum (http://hdfeos.org/forums). 
;
;  If you would like to see an  example of any other NASA HDF/HDF-EOS data 
; product that is not listed in the HDF-EOS Comprehensive Examples page 
; (http://hdfeos.org/zoo), feel free to contact us at eoshelp@hdfgroup.org 
; or post it at the HDF-EOS Forum (http://hdfeos.org/forums).
;
; Acknowledgement:
;
;  We would like to thank Ryan Pavlick for suggestions and bug fixes.
; 
; Tested under: NCL 6.1.0
; Last Updated: 2013-10-28

load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
; Load kmlncl package [1].
load "$NCARG_ROOT/lib/ncarg/nclscripts/kml/kml.ncl"


begin
; To read HDF-EOS5 files, .he5 is required.
; For more information, consult section 4.3.2 of [2].
  file_name = "OMI-Aura_L2-OMNO2_2008m0720t2016-o21357_v003-2008m0721t101450.he5"
  eos_file=addfile(file_name, "r") ; Read file. 

  data_unscaled=eos_file->CloudFraction_ColumnAmountNO2

; Check attributes.
  printVarSummary(data_unscaled)

  scale_factor=data_unscaled@ScaleFactor
  add_offset=data_unscaled@Offset

; Apply scale and offset. The equation is scale*(data-offset).
  data=scale_factor*(data_unscaled-add_offset)

  data@lat2d=eos_file->Latitude_ColumnAmountNO2
  data@lon2d=eos_file->Longitude_ColumnAmountNO2

; Reuse the original attributes for nice plot.
  data@long_name = data_unscaled@Title
  data@units = data_unscaled@Units


; The following resource attributes are slightly modified from the original [3]
; to produce a good KML.
  res=True

; Turn the followings off for correct KML mapping.
; res@gsnMaximize=True 
; res@gsnPaperOrientation = "portrait"  
; res@lbOrientation="vertical" 
; res@tiMainString = file_name


; Use the entire color spectrum.
  res@gsnSpreadColors=True   

; Enable contour fill.
  res@cnFillOn=True   

; This is added for KML. 
  res@mpLandFillColor="transparent" 

 ; Turn off contour line.
  res@cnLinesOn=False  

; Turn off contour line labels.
  res@cnLineLabelsOn =  False 

; This preserves the grid cell shapes from the HDF file. 
; This dramatically improves how your image looks in the Arctic.
; It also speeds up the ncl script.
  res@cnFillMode="CellFill"   
  res@lbLabelAutoStride= True


; Plot type must be PS for kmlncl.
  wks_type = "ps"

; Bigger size means a better resolution on Google Earth.
  wks_type@wkPaperSize="A0"

; Open workstation.
  xwks = gsn_open_wks (wks_type, file_name)
  gsn_define_colormap(xwks,"BlAqGrYeOrReVi200")
  plot=gsn_csm_contour_map(xwks,data,res)

; Open a new KML document.
  kres=True
  kml = kml_open_document(file_name, data@long_name, kres) 

; Position resources for ScreenOverlay element.
  kres@kmlOverlayXY  = True
  kres@kmlOverlayXYx  = 0
  kres@kmlOverlayXYy  = 0
  kres@kmlOverlayXYxunits = "fraction" ;
  kres@kmlOverlayXYyunits = "fraction" ; pixels, or insetPixels
  kres@kmlScreenXY = True
  kres@kmlScreenXYx = 0
  kres@kmlScreenXYy = 0
  kres@kmlScreenXYxunits = "fraction"
  kres@kmlScreenXYyunits = "fraction"
  
; Crop the labelbar from the plot and converts it to a png.
  lbicon = kml_crop_labelbar(kml, "label bar", xwks, plot, kres)
  
; Add labelbar to KML document as a ScreenOverlay element.
  kml = kml_add_screenoverlay (kml, "Label Bar", lbicon, kres )

; Add cropped labelbar filename to list of files to be compressed in KMZ 
; archive.
  kml_add_kmzfile(kml,lbicon)
  
; Get the latlon coordinates describing the corners of the groundoverlay plot.
  LatLonBox = kml_get_vp_latlonbox(plot)

; Get the page coordinates of the actual map from the plot.
  cropbox = kml_get_vp_cropbox(xwks,plot)

; You must "delete(wks)" before cropping map.
  delete(xwks)

; kml_crop_plot() will now set white areas of the map to transparent.
  kres@kmlCropAlphaColor = "white"         

; Crop the groundoverlay map and converts it to png.
  icon = kml_crop_plot(cropbox, kres)      

; Set alpha channel for overlay to make it appear partially transparent in 
; Google Earth.
  kres@kmlColor = "99ffffff"                                                    

; Add cropped map to KML document as a GroundOverlay element.
  kml = kml_add_groundoverlay(kml, "GroundOverlay", icon, LatLonBox, 0, kres)

; Add cropped map filename to list of files to be compressed in KMZ archive.
  kml_add_kmzfile(kml,icon)

  delete(kres@kmlColor)

; Close kml document.
  kml = kml_close_document(kml) 

; Write kml document.
  kml_write(kml)            

; Zip kml file and overlay images together in to kmz.
  kml_make_kmz(kml)         

  delete(plot)
  delete(data)
  delete(res)
  delete(eos_file)
end

; References
;
; [1] http://consideragainthat.org/kmlncl/index.html
; [2] http://hdfeos.org/software/ncl.php
; [3] http://hdfeos.org/zoo/GESDISC/OMI-Aura_L2-OMNO2_2008m0720t2016-o21357_v003-2008m0721t101450_CloudFraction.ncl