代码分享及程序运行时间问题求助。。。

jyfish · 发表于 2017-3-13 17:18:56

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

刚入门NCL，编的程序有点糙，多包涵，实现功能是全球1000多个站点的ASCII（txt）数据插值到格点数据，然后写入到NC文件。听说NCL对循环处理不好，可是从Matlab 过来的，习惯了循环读取文件，然后相应写出到对应文件。实在想不出来更好的解决办法。。。

使用的NCL version 是6.1.2，故 cd_inv_string 和 get_unique_values 还不支持。所以添加了两个子程序实现上述功能，第三个子程序是按照客户需求编的插值算法，可用obj_anal_ic_Wrap 等插值函数替代。

程序运行基本没有问题，提示的warning 不影响运行结果，生成的NC数据没有问题。现在问题来了，第一天的数据用时2分钟，第二天的3分钟，后面依次增加时间，循环递增，时间也在递增，到后面一次要一个小时生成一个nc文件。经过程序运行的时间节点的分析，发现调用子程序运行的时间在不断增加。但我实在不知道解决办法。。。。求大神指点，谢谢！

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
undef("cd_inv_string")
function cd_inv_string(instring, format)
local months, fmonths, century, time, year, month, day, hour, minute,\
second, doy, inchars, strm, c, n, pm_code, monday, units, format, informat
begin
months = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
fmonths = (/"","January","February","March","April","May","June","July","August","September","October","November","December"/)
century = get_res_value_keep(format, "century", 2000)
units = get_res_value_keep(format, "units", "hours since 1800-01-01 00:00:00")
time = new( dimsizes(instring), double)
do n=0, dimsizes(instring)-1
year = 2000
month = 01
day = 01
hour = 00
minute = 00
second = 00
doy = 101
doy@used = False
if(dimsizes(format).eq.dimsizes(instring))
inchars = tochar(format(n))
else
inchars = tochar(format(0))
end if
strm = 0
do c=0, dimsizes(inchars)-1
if(inchars(c) .ne. "%")
strm = strm+ 1
continue
else
c= c+1
pm_code = inchars(c)
if(pm_code .eq. "Y")
year = toint(str_get_cols(instring(n), strm, strm+3 ) )
strm = strm+4
end if
if(pm_code .eq. "y")
year = toint(str_get_cols(instring(n), strm, strm+1 ) ) + century
strm = strm+2
end if
if(pm_code .eq. "N")
month = toint(str_get_cols(instring(n), strm, strm+1 ) )
strm = strm+2
end if
if(pm_code .eq. "D")
day = toint(str_get_cols(instring(n), strm, strm+1 ) )
strm = strm+2
end if
if(pm_code .eq. "H")
hour = toint(str_get_cols(instring(n), strm, strm+1 ) )
strm = strm+2
end if
if(pm_code .eq. "M")
minute = toint(str_get_cols(instring(n), strm, strm+1 ) )
strm = strm+2
end if
if(pm_code .eq. "S")
second = toint(str_get_cols(instring(n), strm, strm+1 ) )
strm = strm+2
end if
if(pm_code .eq. "f")
second = toint(str_get_cols(instring(n), strm, strm+1 ) )*0.6 ;;; *100/60
strm = strm+2
end if
if(pm_code .eq. "J")
doy = toint(str_get_cols(instring(n), strm, strm+2 ) )
doy@used = True
strm = strm+2
end if
end if
end do
if(doy@used) ; put it at the end so doy can come before a year is defined
monday = monthday(year, doy)
month = toint(monday)/100
day = toint( toint(monday) - (month*100) )
end if
time(n) = cd_inv_calendar(year, month, day, hour, minute, second, units,0)
end do
return(time)
end
undef("get_unique_values")
function get_unique_values(var)
local varid, var, levels, i, dummy
begin
var1d = ndtooned(var)
;printVarSummary(var1d)
levels = new(1,typeof(var1d))
do i=0,dimsizes(var1d)-1
if (ismissing(levels(0)) .eq. True) then
if (ismissing(var1d(i)) .ne. True) then
levels(0) = var1d(i)
end if
else
if (any(levels.eq.var1d(i)) .or. ismissing(var1d(i)).eq.True) then
else
dummy = new(dimsizes(levels) + 1, typeof(levels))
dummy(0:dimsizes(levels)-1) = levels
dummy(dimsizes(dummy)-1) = var1d(i)
delete(levels)
levels = dummy
delete(dummy)
end if
end if
end do
return(levels)
end
undef("mean_within_grid")
function mean_within_grid(lon, lat, data, olon, olat)
local lon,lat,data,olon,olat, nlon, nlat, new_data, i, j, \
min_lon, max_lon, min_lat, max_lat, dummy_lat, x_index, \
y_index, dummy_data, index
begin
nlon = dimsizes(olon)
nlat = dimsizes(olat)
new_data = new((/nlat, nlon/), float, -9999)
do i = 0, nlon-1
do j = 0, nlat-1
min_lon = olon(i)-5
max_lon = olon(i)+5
min_lat = olat(j)-5
max_lat = olat(j)+5
y_index = ind((lon.gt.min_lon).and.(lon.le.max_lon))
if any(ismissing(y_index)).eq.False then
dummy_lat = lat(y_index)
x_index = ind((dummy_lat.gt.min_lat).and. \
(dummy_lat.le.max_lat))
if any(ismissing(x_index)).eq.False then
index= y_index(x_index)
dummy_data = data(index)
new_data(j,i) = dim_avg_Wrap(dummy_data)
delete([/index,dummy_data/])
else
new_data(j,i) = -9999
end if
delete([/x_index, dummy_lat/])
else
new_data(j,i) = -9999
end if
delete([/min_lon, max_lon, min_lat, max_lat, \
y_index/])
end do
end do
return(new_data)
end
begin
tstart = "20000101" ;user input
t1 = cd_inv_string(tstart, "%Y%N%D")
tend = "20000110" ;user input
t2 = cd_inv_string(tend, "%Y%N%D")
hours = 0 ;hours options are 0 or 12
common_name = "UPAR_WEA_GLB_MUL_FTM_MERGE_radiosonde_"
common_inpath="/cma/u/Twangzzh/ClimateModel/yjiang/DATAOUT/OUT2000-2016_with_GIS/"
common_outpath = "/cma/u/Twangzzh/ClimateModel/yjiang/NCL/DATAOUT/"
do ntimes = 0,(t2-t1)/24
dummy_time = t1 + 24*ntimes + hours; unit in hours since 1800-01-01 00:00:00
dummy_time@units = "hours since 1800-01-01 00:00:00"
utc_date = cd_calendar(dummy_time, 0)
year = tointeger(utc_date(:,0))
month = tointeger(utc_date(:,1))
day = tointeger(utc_date(:,2))
hour = tointeger(utc_date(:,3))
date_str = sprinti("%0.4i", year) + sprinti("%0.2i", month) + sprinti("%0.2i", day)\
+ "-" + sprinti("%0.2i", hour)
in_path = common_inpath + sprinti("%0.4i", year) + "/" + sprinti("%0.2i", month) + "/"
input = in_path + common_name + date_str + ".txt"
if isfilepresent(input).eq. False then
delete([/dummy_time, utc_date, year, month, day, hour, date_str, in_path, input/])
continue
end if
dummy := asciiread(input, -1, "string")
nrows = numAsciiRow(input)
;printVarSummary(dummy)
;print(nrows)
; if str_is_blank(dummy(nlines-1)).eq.False
; f = dummy
; else
; f = dummy(0:nlines-2)
; end if
f = dummy(0:nrows-1)
;printVarSummary(f)
station_name = str_get_cols(f, 0, 11);
time_str = str_get_cols(f, 12, 24);
format = "%Y %N %D %H"
format@centry = 1900
format@units= "hours since 1900-01-01 00:00:00"
system("date")
time_s = cd_inv_string(time_str, format)
system("date")
time = get_unique_values(time_s)
system("date")
delete(time_s)
time!0 = "time"
time@long_name = "time"
time@units = "hours since 1900-01-01 00:00:00"
time@_FillValue = -9999.
time&time = time
print("time = " + time)
;printVarSummary(time)
lat = stringtofloat(str_get_field(f, 6, " "));
lon = stringtofloat(str_get_field(f, 7, " "));
lon = where(lon.lt.0, lon+360, lon);
pre = stringtofloat(str_get_field(f, 8, " "))*0.01; Pressure in unit pa converting to millibar
Ht = stringtofloat(str_get_field(f, 9, " ")); Height in m
Ht = where(Ht.eq.-9999., -9999., Ht);
Ht@_FillValue = -9999.;
temp = stringtofloat(str_get_field(f, 10, " ")); Temperature in C
temp = where(temp.eq.-9999., -9999., temp*0.1);
temp@_FillValue = -9999.;
TD = stringtofloat(str_get_field(f, 11, " ")); Dew temeprature in C
TD = where(TD.eq.-9999., -9999., TD*0.1);
TD@_FillValue = -9999.;
WD = stringtofloat(str_get_field(f, 12, " ")); Wind direction in degrees north
WD = where((WD.eq.-9999.).or.(WD.eq.999.), -9999., WD);
WD@_FillValue = -9999.;
WS = stringtofloat(str_get_field(f, 13, " ")); Wind speed in m/s
WS = where(WS.eq.-9999., -9999., WS*0.1);
WS@_FillValue = -9999.;
olat = ispan(-85,85,10);
olat!0 = "latitude"
olat@long_name = "latitude"
olat@units ="degrees_north"
olat&latitude = olat
olon = ispan(5,355,10);
olon!0 = "longitude"
olon@long_name = "longitude"
olon@units = "degrees_east"
olon&longitude = olon
PL = (/1,2,3,5,7,10,20,30,50,70,100,125,150,175,\
200,225,250,300,350,400,450,500,550,600,650,700,750,\
775,800,825,850,875,900,925,950,975,1000/); defined pressure levels
PL!0 = "level"
PL@long_name = "pressure level"
PL@units = "millibars"
PL&level = PL
linlog = 1 ;use linear interpolation and avoid out of range extrapolation
; vertical interpolation over pressure levels
station = get_unique_values(station_name)
nstation = dimsizes(station)
new_lat = new(nstation, float)
new_lon = new(nstation, float)
new_temp_data = new((/dimsizes(PL), nstation/), float)
new_TD_data = new((/dimsizes(PL), nstation/), float)
new_WD_data = new((/dimsizes(PL), nstation/), float)
new_WS_data = new((/dimsizes(PL), nstation/), float)
;print(station)
system("date")
print("ntimes="+ntimes+",loop nstation="+(nstation-1))
do ista = 0, nstation-1
;print("ntimes="+ntimes+",ista="+ista)
lat_data = get_unique_values(lat(ind(station_name.eq.station(ista))))
lon_data = get_unique_values(lon(ind(station_name.eq.station(ista))))
pre_data = pre(ind(station_name.eq.station(ista)))
temp_data = temp(ind(station_name.eq.station(ista)))
TD_data = TD(ind(station_name.eq.station(ista)))
WD_data = WD(ind(station_name.eq.station(ista)))
WS_data = WS(ind(station_name.eq.station(ista)))
if dimsizes(get_unique_values(pre_data)).gt.1 then
new_temp = int2p_n(pre_data(::-1), temp_data(::-1), PL, linlog, 0);
new_TD = int2p_n(pre_data(::-1), TD_data(::-1), PL, linlog, 0);
new_WD = int2p_n(pre_data(::-1), WD_data(::-1), PL, linlog, 0);
new_WS = int2p_n(pre_data(::-1), WS_data(::-1), PL, linlog, 0);
else
new_temp = fspan(-9999, -9999, dimsizes(PL));
new_TD = fspan(-9999, -9999, dimsizes(PL));
new_WD = fspan(-9999, -9999, dimsizes(PL));
new_WS = fspan(-9999, -9999, dimsizes(PL));
new_temp@_FillValue = -9999;
new_TD@_FillValue = -9999;
new_WD@_FillValue = -9999;
new_WS@_FillValue = -9999;
end if
;print(PL + " " +new_temp)
new_lat(ista) = lat_data;
new_lon(ista) = lon_data;
new_temp_data(:,ista) = new_temp;
new_TD_data(:,ista) = new_TD;
new_WD_data(:,ista) = new_WD;
new_WS_data(:,ista) = new_WS;
delete([/lat_data, lon_data, pre_data, temp_data, TD_data, WD_data, WS_data,\
new_temp, new_TD, new_WD, new_WS/]);
end do
;print(new_temp_data(:,0)+ " " + new_temp_data(:,1) + " " + new_temp_data(:,2));
;print(new_lon + " " + new_lat + " " + new_temp_data(35,:));
;print(lat_data+ " " +lon_data+ " " +pre_data+ " "+ temp_data);
system("date")
; horiontal interpolation over lat and lon
final_temp = new((/dimsizes(time), dimsizes(PL), dimsizes(olat), dimsizes(olon)/), float)
final_TD = new((/dimsizes(time), dimsizes(PL), dimsizes(olat), dimsizes(olon)/), float)
final_WD = new((/dimsizes(time), dimsizes(PL), dimsizes(olat), dimsizes(olon)/), float)
final_WS = new((/dimsizes(time), dimsizes(PL), dimsizes(olat), dimsizes(olon)/), float)
;rscan = (/100, 70, 40, 10/);
;opt = False;
;final_temp(0,:,:,:) = obj_anal_ic_Wrap(new_lon,new_lat,new_temp_data, olon, olat,rscan, opt)
system("date")
;time dimension here is 1, so s_ntime = 1
s_ntime = dimsizes(time)
do nlevel = 0 , dimsizes(PL)-1
final_temp(s_ntime-1, nlevel,:,:) = mean_within_grid(new_lon, new_lat, new_temp_data(nlevel,:), olon, olat)
final_TD(s_ntime-1, nlevel,:,:) = mean_within_grid(new_lon, new_lat, new_TD_data(nlevel,:), olon, olat)
final_WD(s_ntime-1, nlevel,:,:) = mean_within_grid(new_lon, new_lat, new_WD_data(nlevel,:), olon, olat)
final_WS(s_ntime-1, nlevel,:,:) = mean_within_grid(new_lon, new_lat, new_WS_data(nlevel,:), olon, olat)
end do
system("date")
final_temp!0 = "time"
final_temp!1 = "level"
final_temp!2 = "latitude"
final_temp!3 = "longitude"
final_temp@long_name = "Temperature"
final_temp@units = "Celsius_degrees" ; units?
final_temp&time = time
final_temp&level =PL
final_temp&latitude = olat
final_temp&longitude = olon
;print(final_temp(0,35,:,:))
;printVarSummary(final_temp);
final_TD!0 = "time"
final_TD!1 = "level"
final_TD!2 = "latitude"
final_TD!3 = "longitude"
final_TD@long_name = "Dew Temperature"
final_TD@units = "Celsius_degrees" ; units?
final_TD&time = time
final_TD&level =PL
final_TD&latitude = olat
final_TD&longitude = olon
;printVarSummary(final_TD);
final_WD!0 = "time"
final_WD!1 = "level"
final_WD!2 = "latitude"
final_WD!3 = "longitude"
final_WD@long_name = "Wind Direction"
final_WD@units = "degrees" ; units?
final_WD&time = time
final_WD&level =PL
final_WD&latitude = olat
final_WD&longitude = olon
;printVarSummary(final_WD);
final_WS!0 = "time"
final_WS!1 = "level"
final_WS!2 = "latitude"
final_WS!3 = "longitude"
final_WS@long_name = "Wind Speed"
final_WS@units = "meter per second" ; units?
final_WS&time = time
final_WS&level =PL
final_WS&latitude = olat
final_WS&longitude = olon
;printVarSummary(final_WS);
;;;;;;;;;;;;;;;;;output nc file
out_path = common_outpath + sprinti("%0.4i", year) + "/" + sprinti("%0.2i", month) + "/"
output_name = sprinti("%0.4i", year) + sprinti("%0.2i", month) + sprinti("%0.2i", day)\
+ "." + sprinti("%0.2i", hour) + ".pl"
output = out_path + output_name + ".nc"
if isfilepresent(out_path).eq.False then
cmchar = "mkdir -p " + out_path
system(cmchar)
delete(cmchar)
end if
cnchar = "rm -rf " + output
system(cnchar)
fout = addfile(output, "c")
setfileoption(fout, "DefineMode", True)
fileAtt = True
fileAtt@title = "Test NetCDF File"
fileAtt@creation_data = systemfunc("date")
fileattdef(fout,fileAtt)
dimNames = (/"longitude", "latitude", "level", "time"/)
dimSizes = (/dimsizes(olon), dimsizes(olat), dimsizes(PL), 1/)
dimUnlim = (/False, False, False, False/)
filedimdef(fout, dimNames, dimSizes, dimUnlim)
filevardef(fout, "longitude", typeof(olon), getvardims(olon))
filevardef(fout, "latitude", typeof(olat), getvardims(olat))
filevardef(fout, "level", typeof(PL), getvardims(PL))
filevardef(fout, "time", typeof(time), getvardims(time))
filevardef(fout, "t", typeof(final_temp), getvardims(final_temp))
filevardef(fout, "td", typeof(final_TD), getvardims(final_TD))
filevardef(fout, "wd", typeof(final_WD), getvardims(final_WD))
filevardef(fout, "ws", typeof(final_WS), getvardims(final_WS))
filevarattdef(fout, "longitude", olon)
filevarattdef(fout, "latitude", olat)
filevarattdef(fout, "level", PL)
filevarattdef(fout, "time", time)
filevarattdef(fout, "t", final_temp)
filevarattdef(fout, "td", final_TD)
filevarattdef(fout, "wd", final_WD)
filevarattdef(fout, "ws", final_WS)
fout->longitude = (/olon/)
fout->latitude = (/olat/)
fout->level = (/PL/)
fout->time = (/time/)
fout->t = (/final_temp/)
fout->td = (/final_TD/)
fout->wd = (/final_WD/)
fout->ws = (/final_WS/)
; system("rm -f test.png")
; wks = gsn_open_wks("png", "test")
; gsn_define_colormap(wks, "BlAqGrYeOrRe")
; gsn_define_colormap(wks, "rainbow")
; res = True
; res@gsnAddCyclic = False
; res@cnFillOn = True
; res@cnLinesOn = False
; res@cnLevelSpacingF = 0.5
; res@gsnSpreadColors = True
; res@lbAutoLabelStride = True
; plot = gsn_csm_contour_map(wks, final_temp(0,35,:,:), res)
delete([/dummy_time, utc_date, year, month, day, hour, date_str, in_path, input, f, dummy, nrows, station_name, \
time_str, format, time, lat, lon, pre, Ht, temp, TD, WD, WS, olat, olon, PL, linlog, station, nstation, new_lat, \
new_lon, new_temp_data, new_TD_data, new_WD_data, new_WS_data, final_temp, final_TD, final_WD, final_WS, \
out_path, output_name, output, cnchar, fileAtt, dimNames, dimSizes, dimUnlim, fout/])
system("date")
end do
end

复制代码

ayzqs · 发表于 2017-3-13 20:45:17

多谢楼主的分享精神！

{:5_213:}{:5_213:}

		自动登录	找回密码
密码			立即注册

[其他] 代码分享及程序运行时间问题求助。。。

登录后查看更多精彩内容~

相关帖子