You are here: Home / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / DKRZ NCL EUMETSAT IASI regrid to gaussian T63 grid example

DKRZ NCL EUMETSAT IASI regrid to gaussian T63 grid example

The example script regrid the unstructured EUMETSAT IASI data to a gaussian T63 grid using ESMF.

Example script:

;-----------------------------------------------------------
; remap_IASI_EUMETSAT.ncl:
;
; Remap the IASI EUMETSAT data using ESMF to T63 gaussian grid.
;
; dimensions:
; 	nMeas = 162417 ;
; 	nInstr = 1 ;
; 	nSens = 1 ;
; 	nChannel = 5 ;
; 	nAction = 19 ;
; 	nTest = 15 ;
; 	StringLen = 109 ;
; variables:
; 	float LAT(nMeas) ;
; 		LAT:long_name = "LATITUDE (HIGH ACCURACY)" ;
; 		LAT:units = "DEGREE" ;
; 		LAT:codetable = 5001 ;
; 		LAT:_FillValue_0 = 9.96921e+36f ;
; 	float LON(nMeas) ;
; 		LON:long_name = "LONGITUDE (HIGH ACCURACY)" ;
; 		LON:units = "DEGREE" ;
; 		LON:codetable = 6001 ;
; 		LON:_FillValue_0 = 9.96921e+36f ;
; 	float BT_OBS(nMeas, nChannel) ;
; 		BT_OBS:long_name = "Observed brightness temperature" ;
; 		BT_OBS:units = "K" ;
; 		BT_OBS:codetable = 12063 ;
; 		BT_OBS:_FillValue = 9.96921e+36f ;
;
; 05.02.16  kmf
;-----------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin
  print("")
  
  st = get_cpu_time()                               ;-- for calculating the elapsed CPU time

  diri       = "$HOME/data/EUMETSAT/IASI/"
  fname      = "IASI_004_200907170300-200907170600.nc"
  odir       = "./"
  ofile      =  str_get_field(fname,1,".")+"_remap_to_T63.nc"

;-- output file:  T63 grid definition
  gridx      =  192                                     ;-- number of longitude values
  gridy      =   96                                     ;-- number of latitude values
  gridxinc   =  1.875                                   ;-- longitude increment

;-- gaussian latitude values
  latY = (/ -88.57217 ,    -86.72253 ,    -84.86197 ,    -82.99894 , \
	        -81.13498 ,    -79.27056 ,    -77.40589 ,    -75.54106 ,    -73.67613 , \
	        -71.81113 ,    -69.94608 ,    -68.08099 ,    -66.21587 ,    -64.35073 ,    -62.48557 , \
	        -60.6204  ,    -58.75521 ,    -56.89001 ,    -55.02481 ,    -53.1596  , \
	        -51.29438 ,    -49.42915 ,    -47.56393 ,    -45.69869 ,    -43.83346 , \
	        -41.96822 ,    -40.10298 ,    -38.23774 ,    -36.37249 ,    -34.50724 , \
    	    -32.64199 ,    -30.77674 ,    -28.91149 ,    -27.04624 ,    -25.18099 , \
    	    -23.31573 ,    -21.45048 ,    -19.58522 ,    -17.71996 ,    -15.8547  , \
    	    -13.98945 ,    -12.12419 ,    -10.25893 ,     -8.393669,     -6.528409, \
    	     -4.66315 ,     -2.79789 ,     -0.9326299,     0.9326299,     2.79789 , \
    	      4.66315 ,      6.528409,      8.393669,     10.25893 ,     12.12419 , \
    	     13.98945 ,     15.8547  ,     17.71996 ,     19.58522 ,     21.45048 , \
    	     23.31573 ,     25.18099 ,     27.04624 ,     28.91149 , \
    	     30.77674 ,     32.64199 ,     34.50724 ,     36.37249 ,     38.23774 , \
    	     40.10298 ,     41.96822 ,     43.83346 ,     45.69869 ,     47.56393 , \
    	     49.42915 ,     51.29438 ,     53.1596  ,     55.02481 ,     56.89001 , \
    	     58.75521 ,     60.6204  ,     62.48557 ,     64.35073 ,     66.21587 , \
    	     68.08099 ,     69.94608 ,     71.81113 ,     73.67613 ,     75.54106 , \
    	     77.40589 ,     79.27056 ,     81.13498 , \
    	     82.99894 ,     84.86197 ,     86.72253 ,     88.57217 /)

;-- saver way to define lonX (lonX: 192 grid points, starting at -180.0 equally spaced 1.875)
  lonX = ispan(0,(gridx-1),1)*gridxinc
  lonX = where(lonX .lt. 180.,lonX,lonX-360.)
  qsort(lonX)                                             ;-- sort values to get -180.0-178.125

;-- create coordinate and named coordinate dimensions
  latY!0     = "lat"
  latY&lat   =  latY
  latY@units = "degrees_north"
  lonX!0     = "lon"
  lonX&lon   =  lonX
  lonX@units = "degrees_east"
  
;-- open file
  if (isfilepresent(diri+"/"+fname)) then  
     f = addfile(diri+"/"+fname,"r")
  else
     print("++++  File not found !")
  end if
     
;-- assign data (LON, LAT and BT_OBS) from source netCDF file
  bt1    =  f->BT_OBS(:,2)
  lat1d  =  f->LAT(:)
  lon1d  =  f->LON(:)

  if(any(ismissing(bt1))) then
    print("+++ data contains missing values ...")
  end if

;-- print date to logfile to retrieve run time
  date = systemfunc("date")
  print("--> ESMF regridding start - "+fname+"  "+date)

;-- ESMF resource settings
  reopt                 =  True                           ;-- ESMF resource settings

  reopt@ForceOverwrite  =  True
  reopt@WgtFileName     =  odir+"bt_to_T63.nc"            ;-- weight file
  reopt@NoPETLog        =  True                           ;-- don't generate the "PET0.RegridWeightGen.log" file
  reopt@RemoveSrcFile   =  True                           ;-- remove Src SCRIP grid destination files
  reopt@RemoveDstFile   =  True                           ;-- remove Dst SCRIP grid destination files
    
  reopt@SrcFileName     =  odir+"bt2d_ESMF.nc" ;-- output file
  reopt@DstTitle        = "Swath T63 resolution"          ;-- destination file title
  reopt@DstGridLat      =  latY                           ;-- destination lat grid values
  reopt@DstGridLon      =  lonX                           ;-- destination lon grid values
  reopt@DstRegional     =  True

;  reopt@Debug           =  True
  reopt@SrcRegional     =  True  
  reopt@DstFileName     =  "Swath_T63_SCRIP.nc" ;-- output SCRIP file
  
;-- open a workstation
  wks = gsn_open_wks("png",odir+"plot_IASI_regridded_contour_bilinear_large")

  res                           =  True
  res@gsnAddCyclic              =  False
  res@gsnMaximize               =  True
    
  res@cnFillOn                  =  True  
  res@cnFillMode                = "RasterFill" 
  res@cnLinesOn                 =  False 
  res@cnInfoLabelOn             =  False
  res@cnLineLabelsOn            =  False
  res@cnLevelSelectionMode      = "ManualLevels"
  res@cnMinLevelValF            =  220.
  res@cnMaxLevelValF            =  300.
  res@cnLevelSpacingF           =    2.
  res@cnFillPalette             = "BlAqGrYeOrReVi200"
   
  res@gsnFrame                  = False
  res@gsnDraw                   = False
   
  cnres                         = res
   
  cnres@gsnRightString          = ""
  cnres@gsnLeftString           = ""
  cnres@lbLabelBarOn            = False

  mpres                         = res

  mpres@mpOutlineOn             =  True
  mpres@mpGridAndLimbOn         =  True
  mpres@mpGridLineColor         = "grey30"
  mpres@mpGridLineThicknessF    =  1
  mpres@mpGridLineDashPattern   =  0
  mpres@mpGridAndLimbDrawOrder  = "PostDraw"

  mpres@tiMainString            = "IASI EUMETSAT - ESMF regridded to T63"

;-- regrid the data in rectangular blocks
  lats  = ispan(-90,90,45)          ;-- result rough but the best
  lons  = ispan(-180,180,90)        ;-- result rough but the best
  nlats = dimsizes(lats)
  nlons = dimsizes(lons)

  first_segment = True

  do nlt=0,nlats-2
    do nln=0,nlons-2
      minlat = lats(nlt)
      maxlat = lats(nlt+1)
      minlon = lons(nln)
      maxlon = lons(nln+1)

;-- find all the lat/lon points in this rectangular block
      ii := ind(lat1d.ge.minlat-1.and.lat1d.le.maxlat+1.and.lon1d.ge.minlon-1.and.lon1d.le.maxlon+1)
      print("Range lat="+minlat+":"+maxlat+", lon="+minlon+":"+maxlon)
      if(any(ismissing(ii)).or.dimsizes(ii).le.3) then
         print("No lat/lon values found in this range.")
         continue
      else
         print(dimsizes(ii) + " values found in this range.")
      end if

;-- ESMF regridding
      reopt@SrcGridLat := lat1d(ii)                          ;-- source lat grid values
      reopt@SrcGridLon := lon1d(ii)                          ;-- source lon grid values
      bt1_regrid        = ESMF_regrid(bt1(ii), reopt)        ;-- take some time

;---Create plots of each subset of data.
      if(first_segment) then
         plot_regrid = gsn_csm_contour_map(wks,bt1_regrid,mpres)
         first_segment  = False
      else
         overlay_regrid = gsn_csm_contour(wks,bt1_regrid,cnres)
         overlay(plot_regrid,overlay_regrid)
      end if
    end do
  end do

  draw(plot_regrid)
  frame(wks)

;-- print date to logfile to retrieve run time
  et = get_cpu_time()                               ;-- for calculating the elapsed CPU time
  cputime = et-st
  print("----> CPU time :  "+cputime+" s")


end

Document Actions