[ncl-talk] About regridding from
Tao Lu
hakufu.asano at gmail.com
Wed Jan 18 04:47:13 MST 2017
Hello
Thank you for your reply.
I have succeeded to do objective analysis according to Barry's code.
The problem relied on that obj_anal_ic function does not return meta data.
So, I could not get the plot graph correctly.
Instead of obj_anal_ic, I took the function obj_anal_ic_Wrap.
I attached the edited script.
Thank you for your kind help.
Tao
On Wed, Jan 18, 2017 at 5:56 AM, Barry Lynn <barry.h.lynn at gmail.com> wrote:
> Hi:
> I think it is easiest to grid the radar data to the WRF grid, as I did in
> my original program.
>
> This doesn't use regrid programs, but objective analysis.
>
> If the wrf grid is coarser than the radar data, the data will be smoothed.
>
> Using NCL requires perseverance. Sometimes it is easier to start with a
> basic example.
>
> If you make data lat, Lon, value, the program I sent will read it.
>
> You then need to set/list the files, and the name of the WRF file,
>
> Sorry I can't help more at the moment.
>
>
>
> On Tuesday, January 17, 2017, Dennis Shea <shea at ucar.edu> wrote:
>
>> [1]
>> A general rule for regridding is to go from high-resolution to
>> low-resolution. *If* the variable is 'smoothly varying' virtually any
>> interpolation method (eg: bilinear) will do. If the data is fractal (eg:
>> hourly/3hrly/daily precipitation), then other methods may be more
>> appropriate.
>>
>> Interpolating from low-resolution to high-resolution, does *not* provide
>> any additional information. However, depending upon your objective and use,
>> this may be what you need to do.
>>
>> [2]
>> Look at ESMF Regridding: http://www.ncl.ucar.edu/Applications/ESMF.shtml
>>
>> Look at Example 21
>>
>> [a] Take the radar data and make one-dimensional
>> lat_radar/lon_radar/value_radar arrays. These are your input ('Source').
>> Create a weight file (say)
>>
>> Opt at WgtFileName= "WGT.Radar_to_WRF.nc" ; default is "weights_file.nc"
>>
>> [b] The target grid is defined by the WRF: XLAT, XLONG variables. There
>> are numerous examples on the ESMF page to illustrate how to do this.
>>
>> [c] Run the script
>>
>> [d] The weight file can be used subsequently for other regriddings.
>>
>> =============
>> NOTE:
>> It is your responsibility to look at and study the ESMF URL. Read it
>> carefully before proceeding. Responding to ncl-talk questions, especially,
>> ESMF can be quite time consuming and we have our own jobs to do.
>>
>> Good Luck
>>
>> On Tue, Jan 17, 2017 at 10:41 AM, Barry Lynn <barry.h.lynn at gmail.com>
>> wrote:
>>
>>> Hi Tao:
>>>
>>> I apologize, but I can't answer right now as I am traveling.
>>>
>>> However, in order to provide ncl-talkers with enough information/data,
>>> you can write out the data grid (lat/lon) you are trying to plot on from
>>> WRF, and upload one time step of it (or two). You should also upload the
>>> radar data.
>>>
>>> ;---Write wrf variables lat/lon to file for ncltalk.nc
>>>
>>>
>>> system("rm -f file_for_ncltalk.nc")
>>> fout = addfile("file_for_ncltalk.nc","c")
>>> fout->xlat = xlat
>>> fout->xlong = xlong
>>>
>>> Barry
>>>
>>> On Tue, Jan 17, 2017 at 7:14 PM, Tao Lu <hakufu.asano at gmail.com> wrote:
>>>
>>>> Hello Barry:
>>>>
>>>> Thank you very much for your reply.
>>>> I am very appreciated that you shared your codes with us.
>>>>
>>>> Some editions have been done in my situation.
>>>> However, I could not get contour after grid =
>>>> obj_anal_ic(lon_o,lat_o,obs_o,xlon_m,xlat_m,rscan,False).
>>>> I attached the edited code and run information, please help me check
>>>> that.
>>>>
>>>> I also have some questions about the method for regridding:
>>>> 1. There is a item named "rightmost dimension or leftmost dimension" in
>>>> http://www.ncl.ucar.edu/Document/Functions/Built-in/obj_anal
>>>> _ic.shtml
>>>> I googled it but have no ideal.
>>>>
>>>> 2. Also according to the webpage in 1. obs_o at _FillValue will be
>>>> ignored.
>>>> Does it mean that missing value won't be accounted for when
>>>> interpolation?
>>>>
>>>> 3. "cband radar observations: min=0 max=220" before regridding
>>>> "var_zoom min=0 max=23.5354" after regridding
>>>> After regridding the radar density is much smaller than original
>>>> maximum.
>>>> Is it normal? Maybe as you said rscan is important.
>>>>
>>>> 4. When you use
>>>> grid = obj_anal_ic(lon_o,lat_o,obs_o,xlon_m,xlat_m,rscan,False)
>>>>
>>>> You calculate this
>>>> ;---Get slice of lat/lon data at halfway points.
>>>> nlat = dims2d(0)
>>>> nlon = dims2d(1)
>>>> xlat_m = xlat2d_m(:,nlon/2)
>>>> xlon_m = xlon2d_m(nlat/2,:)
>>>>
>>>> Is it reasonable to calculate lat/lon at a quarter points?
>>>>
>>>> 5. If use halfway points to regrid radar data,
>>>> [image: Inline image 1]
>>>> We will get a n(xlat_m) x n(xlon_m) regridded array.
>>>> So, what will the array's lat/lon will be?
>>>>
>>>> 6. This method is to interpolate radar data to wrf grid.
>>>> Is the mean average method more suitable for radar data? Because
>>>> radar have
>>>> so many data points unlike rain gauge data?
>>>>
>>>>
>>>> My ascii data converted from radar is too large, if you need it I will
>>>> send it to you.
>>>> The data is something like below:
>>>> # crain_composite_data:rain_density:lat:lon:ndata=2406400:
>>>> # 2015-09-09_00:00:00
>>>> -9999 47.32916666666667 137.01874999999998
>>>> -9999 47.32916666666667 137.03125
>>>> -9999 47.32916666666667 137.04375
>>>> -9999 47.32916666666667 137.05625
>>>> -9999 47.32916666666667 137.06875
>>>> -9999 47.32916666666667 137.08124999999998
>>>>
>>>>
>>>> Waiting for your reply.
>>>> Thank you very much.
>>>>
>>>> Tao
>>>>
>>>>
>>>> On Tue, Jan 17, 2017 at 4:05 AM, Barry Lynn <barry.h.lynn at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi Tao:
>>>>>
>>>>> Here, I gridded the radar data to the WRF grid.
>>>>>
>>>>> You will need to adapt what I did to your situation, but the most
>>>>> important step is to decide on the degree of smoothing (which is set by
>>>>> rscan).
>>>>>
>>>>> grid = obj_anal_ic(lon_o,lat_o,obs_o,xlon_m,xlat_m,rscan,False)
>>>>>
>>>>> Please write if you have questions.
>>>>>
>>>>> Barry
>>>>>
>>>>> On Mon, Jan 16, 2017 at 9:11 AM, Tao Lu <hakufu.asano at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Dear ncl-talk
>>>>>>
>>>>>> Hello.
>>>>>>
>>>>>> I am new to regridding, if possible please give me some help and
>>>>>> advice.
>>>>>>
>>>>>> What I am doing?: Compare two kind of grid data
>>>>>> *I have two kind of data.*
>>>>>> 1. radar data:
>>>>>> resolution: 0.5 minute x 0.75 minute, as seen not regular
>>>>>> grids
>>>>>> file format: ascii
>>>>>> [image: Inline image 1]
>>>>>> 2. output data from WRF:
>>>>>> dx = 30000,
>>>>>> dy = 30000,
>>>>>> map_proj = 'lambert',
>>>>>> ref_lat = 36.00,
>>>>>> ref_lon = 140.00,
>>>>>> truelat1 = 30.0,
>>>>>> truelat2 = 60.0,
>>>>>> stand_lon = 140.0,
>>>>>> file fomat: netcdf
>>>>>> [image: Inline image 3]
>>>>>>
>>>>>> *I have to regrid the two data to each other or to new grid.*
>>>>>> 1. Regrid wrf data grid to radar data grid?
>>>>>> 2. Regrid radar data to wrf data?
>>>>>> 3. Regrid both to another data?
>>>>>> Which is easiest?
>>>>>>
>>>>>> I am really new to this, please give me some guide line, I will
>>>>>> further operation on my work.
>>>>>>
>>>>>> Thank you very much for your kind help.
>>>>>>
>>>>>> Tao
>>>>>>
>>>>>>
>>>>>> --
>>>>>> ******************************************************
>>>>>> Tao Lu (M. Eng)
>>>>>> Laboratory of River Engineering and Hydrology,
>>>>>> Dept. of Civil and Environmental Engineering,
>>>>>> Graduate School of Science and Engineering,
>>>>>> Chuo University
>>>>>>
>>>>>> 1-13-27,Kasuga,Bunkyo-ku,Tokyo
>>>>>> 112-8551, Japan
>>>>>> TEL: 03-3817-1805; Phone: 070-2188-7509
>>>>>> Email1: hakufu.asano at gmail.com
>>>>>> Email2: lutao at civil.chuo-u.ac.jp
>>>>>> *****************************************************
>>>>>>
>>>>>> _______________________________________________
>>>>>> ncl-talk mailing list
>>>>>> ncl-talk at ucar.edu
>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Barry H. Lynn, Ph.D
>>>>> Senior Lecturer,
>>>>> The Institute of the Earth Science,
>>>>> The Hebrew University of Jerusalem,
>>>>> Givat Ram, Jerusalem 91904, Israel
>>>>> Tel: 972 547 231 170
>>>>> Fax: (972)-25662581
>>>>>
>>>>> C.E.O, Weather It Is, LTD
>>>>> Weather and Climate Focus
>>>>> http://weather-it-is.com
>>>>> Jerusalem, Israel
>>>>> Local: 02 930 9525
>>>>> Cell: 054 7 231 170
>>>>> Int-IS: x972 2 930 9525
>>>>> US 914 432 3108 <(914)%20432-3108>
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> ******************************************************
>>>> 盧 涛 (ル タオ) 修士課程2年
>>>>
>>>> 〒112-8551 東京都文京区春日1-13-27
>>>>
>>>> 中央大学理工学研究科都市環境学専攻
>>>>
>>>> 河川・水文研究室（山田正教授）
>>>>
>>>>
>>>> TEL: 03-3817-3406; Phone: 070-2188-7509
>>>> Email1: hakufu.asano at gmail.com
>>>> Email2: lutao at civil.chuo-u.ac.jp
>>>> *******************************************************
>>>> ******************************************************
>>>> Tao Lu (M. Eng)
>>>> Laboratory of River Engineering and Hydrology,
>>>> Dept. of Civil and Environmental Engineering,
>>>> Graduate School of Science and Engineering,
>>>> Chuo University
>>>>
>>>> 1-13-27,Kasuga,Bunkyo-ku,Tokyo
>>>> 112-8551, Japan
>>>> TEL: 03-3817-1805; Phone: 070-2188-7509
>>>> Email1: hakufu.asano at gmail.com
>>>> Email2: lutao at civil.chuo-u.ac.jp
>>>> *****************************************************
>>>>
>>>
>>>
>>>
>>> --
>>> Barry H. Lynn, Ph.D
>>> Senior Lecturer,
>>> The Institute of the Earth Science,
>>> The Hebrew University of Jerusalem,
>>> Givat Ram, Jerusalem 91904, Israel
>>> Tel: 972 547 231 170
>>> Fax: (972)-25662581
>>>
>>> C.E.O, Weather It Is, LTD
>>> Weather and Climate Focus
>>> http://weather-it-is.com
>>> Jerusalem, Israel
>>> Local: 02 930 9525
>>> Cell: 054 7 231 170
>>> Int-IS: x972 2 930 9525
>>> US 914 432 3108 <(914)%20432-3108>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>
>
> --
> Barry H. Lynn, Ph.D
> Senior Lecturer,
> The Institute of the Earth Science,
> The Hebrew University of Jerusalem,
> Givat Ram, Jerusalem 91904, Israel
> Tel: 972 547 231 170
> Fax: (972)-25662581
>
> C.E.O, Weather It Is, LTD
> Weather and Climate Focus
> http://weather-it-is.com
> Jerusalem, Israel
> Local: 02 930 9525
> Cell: 054 7 231 170
> Int-IS: x972 2 930 9525
> US 914 432 3108 <(914)%20432-3108>
>
>
--
******************************************************
盧 涛 (ル タオ) 修士課程2年
〒112-8551 東京都文京区春日1-13-27
中央大学理工学研究科都市環境学専攻
河川・水文研究室（山田正教授）
TEL: 03-3817-3406; Phone: 070-2188-7509
Email1: hakufu.asano at gmail.com <mail%3Amet.yamos at gmail.com>
Email2: lutao at civil.chuo-u.ac.jp <mail%3Ayamoto at civil.chuo-u.ac.jp>
*******************************************************
******************************************************
Tao Lu (M. Eng)
Laboratory of River Engineering and Hydrology,
Dept. of Civil and Environmental Engineering,
Graduate School of Science and Engineering,
Chuo University
1-13-27,Kasuga,Bunkyo-ku,Tokyo
112-8551, Japan
TEL: 03-3817-1805; Phone: 070-2188-7509
Email1: hakufu.asano at gmail.com <mail%3Amet.yamos at gmail.com>
Email2: lutao at civil.chuo-u.ac.jp <mail%3Ayamoto at civil.chuo-u.ac.jp>
*****************************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170118/7beb155e/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 4669 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170118/7beb155e/attachment-0003.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 392777 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170118/7beb155e/attachment-0004.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 44666 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170118/7beb155e/attachment-0005.png
-------------- next part --------------
;###############################################################
;
; Program to do objective analysis from radar data to wrf grid
;
;
;##############################################################
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
;---Set bounding box
; bounding box for cover kinu watershed
;lats = (/ 30.00, 40.00 /)
;lons = (/ 132.00, 142.00 /)
lats = (/ 28.00, 44.00 /)
lons = (/ 128.00, 144.00 /)
;---We generate plots, but what kind do we prefer?
type = "png"
type at wkWidth = 1500 ; Increase size for a slightly
type at wkHeight = 1500
slice_loc = "1-2" ; 1/2 of bounding box
rscan_plan = "0.25-0.1" ; latitude degrees
box = "_bigbox_"
output_name = "after_obs" + box + "_" + slice_loc + "_" + rscan_plan
wks = gsn_open_wks(type,output_name)
;---Set some basic resources
res = True
res at gsnMaximize = True
res at gsnAddCyclic = False ; Don't add cyclic point
; map resources
res at mpFillOn = True
res at mpOutlineOn = True
res at mpDataBaseVersion = "HighRes"
;res at mpDataBaseVersion = "MediumRes"
; set projection coordinated with wrf
; map sources
;res at tfDoNDCOverlay = True
res at mpLimitMode = "LatLon"
res at mpMinLatF = 28.0
res at mpMaxLatF = 46.0
res at mpMinLonF = 126.0
res at mpMaxLonF = 146.0
res at mpProjection = "LambertConformal"
res at mpLambertParallel1F = 30.
res at mpLambertParallel2F = 60.
res at mpLambertMeridianF = 136.
res at pmTickMarkDisplayMode = "Always"
res at tmXTLabelsOn = False
res at tmYRLabelsOn = False
; label bar sources
res at lbAutoManage = False
res at pmLabelBarHeightF = 0.12
; contour resources
res at cnFillOn = True
res at cnLinesOn = False
res at cnLevelSelectionMode = "ExplicitLevels"
; color of lab cband
res at cnLevels = (/0.1, 1, 5, 10, 20, 50, 100, 210/)
res at cnFillColors = (/"white","cyan", "blue","green","yellow","darkorange", \
"coral1","red","blueviolet"/)
;---Get xlat,xlon from model for interpolation
filepath = "./wrfout_d01_2015-09-08_12:00:00"
a = addfile(filepath,"r")
time = 0
xlat2d_m = wrf_user_getvar(a,"XLAT",time)
xlon2d_m = wrf_user_getvar(a,"XLONG",time)
dims2d = dimsizes(xlat2d_m)
;printVarSummary(xlat2d_m)
;printVarSummary(xlon2d_m)
;printVarSummary(xland)
;print(dims2d)
;print (xlat2d_m)
;---Set zoomin indices
; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y)
; subtract one since we want to use it as an index in NCL
loc = wrf_user_ll_to_ij(a, lons, lats, True)
loc = loc - 1
x_start = loc(0,0)
x_end = loc(0,1)
y_start = loc(1,0)
y_end = loc(1,1)
print ("---Bounding Box Information---")
print ("x_start = " + x_start + " x_end = " + x_end)
print ("y_start = " + y_start + " y_end = " + y_end)
print ("------------------------------")
print ("")
;---Loop of time
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
do itime = 0, 0
if (itime.eq.0)then
asc_file = "./crain_composite_2015.09.09.00.00_UTC"
end if
;---Read ascii data
stime = get_cpu_time()
data = asciiread(asc_file, -1, "float")
header = readAsciiHead(asc_file,2)
date = str_get_field(header(1),2," ")
; read data
radar = data(7::3)
lat_obs = data(8::3)
lon_obs = data(9::3)
; these long_name and units attributes are not necessary. I
; did this so I got better output from print statements below.
radar at long_name = "cband radar observations"
radar at units = "mm/h"
lat_obs at long_name = "observed latitudes"
lat_obs at units = "degrees_north"
lon_obs at long_name = "observed longitudes"
lon_obs at units = "degrees_east"
radar at _FillValue = -999
etime = get_cpu_time()
print("Reading the radar data took " + (etime-stime) + " seconds.")
print ("")
;printVarSummary(radar)
;printVarSummary(lat_obs)
;printVarSummary(lon_obs)
printMinMax(radar,0)
printMinMax(lat_obs,0)
printMinMax(lon_obs,0)
print ("")
;---Get slice of lat/lon data at halfway points.
nlat = dims2d(0)
nlon = dims2d(1)
xlat_m = xlat2d_m(:,nlon/2)
xlon_m = xlon2d_m(nlat/2,:)
;printVarSummary(xlat_m)
;printVarSummary(xlon_m)
;printMinMax(xlat_m,0)
;printMinMax(xlon_m,0)
;---Objective analysis function
; rscan = (/1,0.5,0.25/)
; rscan = (/0.5,0.25,0.225/)
; rscan = (/0.25,0.125,0.05/)
rscan = (/0.25,0.1/)
; rscan = (/1.0/)
;
; obj_anal_ic_Wrap: iterative improvement objective analysis and returns meta data
;
stime = get_cpu_time()
grid = obj_anal_ic_Wrap(lon_obs,lat_obs,radar,xlon_m,xlat_m,rscan,False)
etime = get_cpu_time()
print("objective analysis took " + (etime-stime) + " seconds.")
print ("")
lon_2d_zoom = xlon2d_m(y_start:y_end,x_start:x_end)
lat_2d_zoom = xlat2d_m(y_start:y_end,x_start:x_end)
var_zoom = grid(y_start:y_end,x_start:x_end)
var_zoom at _FillValue = -999
;printVarSummary(grid)
;printMinMax(grid,0)
;printVarSummary(lon_2d_zoom)
;printMinMax (lon_2d_zoom, 0)
;printVarSummary(lat_2d_zoom)
;printMinMax (lat_2d_zoom, 0)
;printVarSummary(var_zoom)
;printMinMax(var_zoom,0)
;print (var_zoom&south_north) ; extract lat and lon
;
;------------------ function without meta data use ----------------------
;------------------ lon_2d_zoom, lat_2d_zoom to plot --------------------
;
;grid = obj_anal_ic(lon_obs,lat_obs,radar,xlon_m,xlat_m,rscan,False)
;var_zoom = grid(y_start:y_end,x_start:x_end)
;var_zoom at _FillValue = -999
; change multi dimension array to one dimension array
lon_1d = ndtooned(lon_2d_zoom)
lat_1d = ndtooned(lat_2d_zoom)
var_zoom_1d = ndtooned(var_zoom)
;printVarSummary(var_zoom_1d)
; by setting the ScalarField resources sfXArray and sfYArray to
; the 1-d or n-d arrays lon and lat,
; to explicitly define the range for the X and Y axes.
;res at sfXArray = lon_2d_zoom ; Associates lat/lon locations
;res at sfYArray = lat_2d_zoom ; with data values.
;--------------------------------------------------------------------------
;---Draw raster contours
stime = get_cpu_time()
;set titles, size and resources
mainstring = date + " (UTC)"
leftstring = "Cband rain after regridding"
rightstring = "dt = 5 minutes"
size = 0.02
res at tiMainString = mainstring
res at gsnLeftString = leftstring ; add the gsn titles
res at gsnRightString = rightstring
res at tiMainFontHeightF = 0.030
res at gsnLeftStringFontHeightF = size ; instead of using txFontHeightF or gsnStringFontHeightF
;res at gsnCenterStringFontHeightF = size ; to set the gsnLeft/Center/RightString font heights,
res at gsnRightStringFontHeightF = size ; individually set each string's font height.
res at cnFillMode = "RasterFill"
; make plot
contour_pre = gsn_csm_contour_map(wks,var_zoom,res)
etime = get_cpu_time()
print("raster contours took " + (etime-stime) + " seconds.")
print ("")
;---Output grid data after objective analysis
;stime = get_cpu_time()
write_table("after_obs" + output_name + ".txt","w",[/lat_1d,lon_1d,var_zoom_1d/],"%8.3f %8.3f %8.3f")
;etime = get_cpu_time()
;print("write grid data after objective took " + (etime-stime) + " seconds.")
;---End loop of time
end do
end
More information about the ncl-talk
mailing list