DKRZ NCL EUMETSAT IASI regrid to gaussian T63 grid example#
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
Result:
