[ncl-talk] ESMF regrid
Karin Meier-Fleischer
meier-fleischer at dkrz.de
Thu Jun 23 23:42:22 MDT 2016
Hi Laura,
I've attached the regridding script for the EUMETSAT IASI data which is
similar to your data if I remember it correct.
The trick is doing the regridding in parts for a better result.
Hope this will help,
Karin
Am 21.06.16 um 22:30 schrieb Laura Fowler:
> Hi Karin:
>
> Thanks for replying to the e-mail that I posted last week-end. I am
> familiar with the examples that you suggested since I am a MPAS
> developer. I think that the issue is that since I am trying to regrid
> a satellite orbit, the source and destination data do have a lot of
> missing grid-points. That is the reason I was looking at example
> ESMF_regrid_8.ncl but I cannot make it work so far. Let me know if you
> have other suggestions, of course.
>
> Cheers,
> Laura
>
>
>
>
> On Mon, Jun 20, 2016 at 1:42 AM, Karin Meier-Fleischer
> <meier-fleischer at dkrz.de> wrote:
>> Hi Laura,
>>
>> you can take the examples to regrid unstructured data to a rectilinear
>> grid like the examples ESMF_regrid_10.ncl, ESMF_regrid_14.ncl and
>> ESMF_regrid_21.ncl.
>>
>> Bye,
>> Karin
>>
>> Am 18.06.16 um 21:58 schrieb Laura Fowler:
>>> Hello:
>>>
>>> What is the best example that is available from the page
>>> http://www.ncl.ucar.edu/Applications/ESMF.shtml to regrid a satellite
>>> orbit to a rectangular grid using ESMF_regrid? I tried to use
>>> ESMF_regrid_8.ncl but unsuccessfully. Is it even possible to do this
>>> using ESMF_regrid?
>>>
>>> I am trying to regrid some orbital data from the CERES satellite into
>>> a rectangular 0.05 deg grid. I tried to include the options
>>> SrcRegional, SrcGridMask, and SrcGridType, but this did not work very
>>> well. I ran across somebody wanting to do the same with TRMM data but
>>> there was not any attached script in his or her post.
>>>
>>> I am trying to regrid data that I plotted on yellowstone (on
>>> yellowstone, see the pdf
>>> plot_CERES_SSF_XTRK-MODIS_Edition4A_2014011000-2014011023.pdf in
>>> /glade/p/work/laura/DATA/CERESdata/CERESorbits.2014-01-10_to_2014-01-14).
>>> My tentative script is in the script regrid.ncl.
>>>
>>> As always, appreciate your help. Thank you,
>>> Laura
>>>
>>>
>> --
>> Dipl. Geophys. Karin Meier-Fleischer
>> Visualization, NCL
>> Application Support
>>
>> Deutsches Klimarechenzentrum GmbH (DKRZ)
>> Bundesstrasse 45a - D20146 Hamburg - Germany
>>
>> Phone: +49 (0)40 460094 126
>> Fax: +49 (0)40 460094 270
>> E-Mail: meier-fleischer at dkrz.de
>> URL: www.dkrz.de
>>
>> Geschäftsführer: Prof. Dr. Thomas Ludwig
>> Sitz der Gesellschaft: Hamburg
>> Amtsgericht Hamburg HRB 39784
>>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
;-----------------------------------------------------------
; 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 at units = "degrees_north"
lonX!0 = "lon"
lonX&lon = lonX
lonX at 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 at ForceOverwrite = True
reopt at WgtFileName = odir+"bt_to_T63.nc" ;-- weight file
reopt at NoPETLog = True ;-- don't generate the "PET0.RegridWeightGen.log" file
reopt at RemoveSrcFile = True ;-- remove Src SCRIP grid destination files
reopt at RemoveDstFile = True ;-- remove Dst SCRIP grid destination files
reopt at SrcFileName = odir+"bt2d_ESMF.nc" ;-- output file
reopt at DstTitle = "Swath T63 resolution" ;-- destination file title
reopt at DstGridLat = latY ;-- destination lat grid values
reopt at DstGridLon = lonX ;-- destination lon grid values
reopt at DstRegional = True
; reopt at Debug = True
reopt at SrcRegional = True
reopt at 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 at gsnAddCyclic = False
res at gsnMaximize = True
res at cnFillOn = True
res at cnFillMode = "RasterFill"
res at cnLinesOn = False
res at cnInfoLabelOn = False
res at cnLineLabelsOn = False
res at cnLevelSelectionMode = "ManualLevels"
res at cnMinLevelValF = 220.
res at cnMaxLevelValF = 300.
res at cnLevelSpacingF = 2.
res at cnFillPalette = "BlAqGrYeOrReVi200"
res at gsnFrame = False
res at gsnDraw = False
cnres = res
cnres at gsnRightString = ""
cnres at gsnLeftString = ""
cnres at lbLabelBarOn = False
mpres = res
mpres at mpOutlineOn = True
mpres at mpGridAndLimbOn = True
mpres at mpGridLineColor = "grey30"
mpres at mpGridLineThicknessF = 1
mpres at mpGridLineDashPattern = 0
mpres at mpGridAndLimbDrawOrder = "PostDraw"
mpres at 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 at SrcGridLat := lat1d(ii) ;-- source lat grid values
reopt at 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
More information about the ncl-talk
mailing list