[ncl-talk] Missing values after ESMF regridding
Bian Qingyun
bianqy at tea.ac.cn
Thu Mar 15 02:08:08 MDT 2018
Hello,
I'm regridding AMSR-E L3 daily snow following the examples on NCL website http://www.ncl.ucar.edu/Applications/Scripts/ESMF_regrid_13.ncl and hdf website http://hdfeos.org/zoo/index_openNSIDC_Examples.php#AMSR_E.
All is fine,but there are missing values along 90°E and 90°W after regridding while there are no missing values in original data along and near 90°E and 90°W 。
Below is my script. I also attached the AMSR-E data file in case.
----------------------
----------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
;-----read data
srcFileName = "AMSR_E_L3_DailySnow_V09_20100311.hdf" ; Source file
eos_file = addfile(srcFileName + ".he2","r")
lon2d = eos_file->GridLon_Northern_Hemisphere
lat2d = eos_file->GridLat_Northern_Hemisphere
; Process NAN (1e51) fillValue for LAMAZ geolocation.
lon2d=where(lon2d .gt. 1000, -999.0, lon2d)
lat2d=where(lat2d .gt. 1000, -999.0, lat2d)
lon2d at _FillValue = -999.0
lat2d at _FillValue = -999.0
hdf_file = addfile(srcFileName, "r")
; Read the dataset.
hdf_data = hdf_file->SWE_NorthernDaily
; Filter out invalid range values.
; See "Table 2. Pixel Values ofr the SWE Feids" from [2].
hdf_data = where(hdf_data .gt. 240, hdf_data at _FillValue, hdf_data)
; Prepare data for plotting by converting type.
swe2d = tofloat(hdf_data)
swe2d at _FillValue = tofloat(hdf_data at _FillValue)
; Multiply by two according to data spec [2].
swe2d = 2 * swe2d
; You can get the description of data set from the data spec [2].
swe2d at long_name = "Northern Hemisphere 5-day Snow Water Equivalent ~C~ (" + hdf_data at hdf_name + ")"
swe2d at units = "mm"
; Associate data with lat/lon.
swe2d at lon2d = lon2d
swe2d at lat2d = lat2d
;-----ESMF regridding
Opt = True ; Options for regridding
Opt at SrcFileName = "EASE_SCRIP.nc" ; Output files
Opt at DstFileName = "NH_SCRIP.nc"
Opt at WgtFileName = "EASE_2_NH_patch.nc"
Opt at ForceOverwrite = True
Opt at SrcMask2D = where(ismissing(lat2d),0,1)
Opt at SrcGridLat = lat2d
Opt at SrcGridLon = lon2d
Opt at DstGridType = "0.25deg" ; Destination grid
Opt at DstTitle = "Northern Hemisphere 0.25 resolution"
Opt at DstLLCorner = (/ 0.25d, 0.25d/)
Opt at DstURCorner = (/89.75d, 179.75d/)
Opt at InterpMethod = "bilinear" ; Careful! Patch method takes a long time
Opt at Debug = True
swe_regrid = ESMF_regrid(swe2d,Opt) ; Regrid swe
printVarSummary(swe_regrid)
;----------------------------------------------------------------------
; Plotting section
;----------------------------------------------------------------------
wks = gsn_open_wks("x11","ESMF_regrid") ; send graphics to PNG file
res = True ; Plot mods desired.
res at gsnMaximize = True ; Maximize plot
res at gsnDraw = False
res at gsnFrame = False
res at cnFillOn = True ; color plot desired
res at cnFillPalette = "amwg" ; set color map
res at cnLinesOn = False ; turn off contour lines
res at cnLineLabelsOn = False ; turn off contour labels
res at cnFillMode = "RasterFill" ; turn raster on
res at cnLevelSelectionMode= "ExplicitLevels" ; set explicit contour levels
res at cnLevels = (/-300,-250,-200,-150,-100, \
0,1,5,10,25,100,200,300,400/)
res at lbLabelBarOn = False ; turn on in panel
res at trGridType = "TriangularMesh" ; allow missing coordinates
res at gsnPolar = "NH" ; specify the hemisphere
res at mpMinLatF = 35
;---Plot original data.
res at gsnAddCyclic = False
res at sfXArray = lon2d
res at sfYArray = lat2d
res at tiMainString = "Original EASE grid (" + str_join(dimsizes(lat2d),",") + ")"
res at trGridType = "TriangularMesh"
res at cnFillMode = "RasterFill"
plot_orig = gsn_csm_contour_map_polar(wks,swe2d,res)
delete(res at sfXArray)
delete(res at sfYArray)
;---Plot regridded data.
res at gsnAddCyclic = True
dims = tostring(dimsizes(swe_regrid))
res at tiMainString = "Regridded to 0.25 degree grid (" + \
str_join(dims," x ") + ")"
plot_regrid = gsn_csm_contour_map_polar(wks,swe_regrid,res)
;---Compare the plots in a panel
pres = True
pres at gsnMaximize = True
pres at gsnPanelLabelBar = True
pres at lbLabelFontHeightF = 0.01
pres at pmLabelBarWidthF = 0.8
gsn_panel(wks,(/plot_orig,plot_regrid/),(/1,2/),pres)
end
---------------------------
---------------------------
Any suggestions are appreciated! Thank you!
Best,
Qingyun
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180315/58969264/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: AMSR_E_L3_DailySnow_V09_20100311.hdf
Type: application/octet-stream
Size: 2156226 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180315/58969264/attachment-0001.obj>
More information about the ncl-talk
mailing list