[ncl-talk] GCOM-C L2 product visualization

[Tomoko Koyama]

Hello.

I'd like to visualize GCOM-C Level 2 tile data which format is hdf5.

The latitude/longitude to each pixel need to be calculated and I wrote a code for visualization.
I followed page 4-15 of the provided manual, GCOM-C "SHIKISAI" Data Users Handbook.
(https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/GCOM-C_SHIKISAI_Data_Users_Handbook_en.pdf)

Unfortunately, the result doesn't look right: I plotted the "over land aerosol optical thickness data", but appeared data are over ocean.
Is this a projection issue? 

If there's any sample code for GCOM-C data visualization, can anybody share with me?

Thank you very much in advance,
Tomoko

;========================================================================
; read/plot GCOM-C Level 2 product (ARNP)

begin
;************************************************
; read in data
;************************************************
  setfileoption("h5", "FileStructure", "Advanced")

  fn = "/Users/tkoyama/Research/SGLI/L2/GC1SG1_20200105D01D_T0529_L2SG_ARNPK_2000.h5"

  var = "AROT_land"
  ;var = "AROT_ocean"

  pltDir	= "/Users/tkoyama/Research/Images/"
  pltName	= "test_"+var
  pltType	= "X11"

  fi = addfile(fn, "r")

  g1 = fi=>/Geometry_data

  ; Get attribute names of the group
  g1atts= getvaratts(g1)
  ;print(g1atts)

  ;---------------------------------------------------------------------------------------
  ; Derive lat/lon of each grid point
  ;
  ; Four corners of "T0529" tile
  ; lat = (/ 30.0, 30.0, 40.0, 40.0 /)
  ; lon = (/127.0171, 138.5641, 143.5948, 156.6489/)
  ;
  pi 		= 3.1416

  ; Folowing page 4-15 of the provided manual,
  ; GCOM-C "SHIKISAI" Data Users Handbook
  ; https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-C/GCOM-C_SHIKISAI_Data_Users_Handbook_en.pdf

  lin_tile 	= 1200
  col_tile 	= 1200
  vtile		= 5
  htile		= 29
  vtile_num	= 18
  htile_num	= 36

  lat = new((/col_tile, lin_tile/), float)
  lon = new((/col_tile, lin_tile/), float)

  d 	= 180.0/lin_tile/vtile_num
  NL 	= 180.0/d
  NP0	= 2*round(NL,3)

  lin_tot	= ispan(0, lin_tile-1, 1) + vtile * lin_tile
  col_tot	= ispan(0, col_tile-1, 1) + htile * col_tile

  lat_tmp 	= 90.0 - (lin_tot + 0.5) * d	; in degree

  NPi	= NP0*cos(lat_tmp*pi/180.0)

  do ilin = 0, lin_tile-1
    lat(:,ilin) = lat_tmp(ilin)
    do jcol = 0, col_tile-1
       lon(jcol,ilin) = 360.0/NPi(ilin) * (col_tot(jcol) - NP0/2 + 0.5)	; in degre
    end do
  end do

  ;---------------------------------------------------------------------------------------
  time = short2flt_hdf( g1->Obs_time )
  t =  g1->Obs_time
  t_atts = getvaratts(t)
  ;print(t@$t_atts(1)$)			; This shows the attributes of the target variable

  ;---------------------------------------------------------------------------------------
  g2 = fi=>/Image_data
  g2atts = getvaratts(g2)

  x = g2->$var$
  ;print(getvaratts(x))

  ; get AOT
  val = where(x.ne.x at Error_DN, x*x at Slope + x at Offset, -999.9)
  val at Spatial_resolution = x at Spatial_resolution
  val at Spatial_resolution_unit = x at Spatial_resolution_unit
  val at Unit = x at Unit
  val at _FillValue	= -999.9
  ;printMinMax(val, True)

;************************************************
; Plot
;************************************************
  wks = gsn_open_wks(pltType, pltDir+pltName)

  setvalues NhlGetWorkspaceObjectId()
      "wsMaximumSize": 100000000	; need some extra workspace
  end setvalues

  res			= True
  res at gsnAddCyclic	= False		; data not global

  res at gsnMaximize	= True		; make ps/eps/pdf large
  res at gsnPaperOrientation	= "portrait"

  res at cnFillOn		= True		; turn on color fill
  res at cnFillPalette	= "matlab_jet"	; set color map
  res at cnLinesOn		= False		; turn off contour lines
  res at cnFillMode	= "RasterFill"	; Raster Mode

  res at cnLevelSelectionMode      = "ManualLevels"
  res at cnMinLevelValF  = 0               ; set the minimum contour level
  res at cnMaxLevelValF  = 2               ; set the maximum contour level
  res at cnLevelSpacingF = 0.2		; set the interval between contours

  res at lbOrientation	= "vertical"	; vertical label barb's
  res at lbLabelFontHeightF= 0.012		; change font size [make smaller]
  res at pmLabelBarWidthF	= 0.1		; make thinner

  res at tiMainString	= var

  res at sfYArray		= lat
  res at sfXArray		= lon

  res at mpMinLatF		= 28.0
  res at mpMaxLatF		= 42.0
  res at mpMinLonF		= 125.0
  res at mpMaxLonF		= 158.0

  ;plot	= gsn_csm_contour(wks,val,res)
  plot	= gsn_csm_contour_map(wks,val,res)

end