Sie sind hier: Startseite / Services / Data Analysis and Visualization / Visualization / Software / NCL / examples / source_code / NCL script to regrid EUMETSAT IASI data to T63 gaussian grid
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

NCL script to regrid EUMETSAT IASI data to T63 gaussian grid

NCL script to regrid EUMETSAT IASI data to T63 gaussian grid.

Plain Text icon NCL_remap_IASI_EUMETSAT_T63_triplets_subsets.ncl — Plain Text, 7 KB (8166 bytes)

Dateiinhalt

;-----------------------------------------------------------
; 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

Artikelaktionen