Dear Eric and ncl-users,
each grid point has a value. For a 3-d time series (nlev x nlat x mlon) the
(0,0,0) grid point has the value 48.955, (0,1,0) 54.956, e.t.c.
Using the onedtond function both (0,0,0) and (0,1,0) give me the same value,
48.955.
I wonder if there is a way to overcome this problem..
Thanks a lot..
Ioannis
From: Rick Brownrigg [mailto:brownrig@ucar.edu]
Sent: Thursday, July 25, 2013 6:45 PM
To: Ioannis Koletsis
Cc: ncl-talk@ucar.edu
Subject: Re: Onedtond problem
Hi,
As near as I can tell, "windprop" is a scalar value, not an array. So
onedtond would simply fill whatever sized array with that one value each
time through the loops. I'm a bit confused as to what you mean by "The
desired results would be (0,0,0) 1st value, (0,1,0) 2nd etc" ?
Rick
On Jul 25, 2013, at 7:17 AM, "Ioannis Koletsis" <koletsis@noa.gr> wrote:
Dear ncl users,
I calculate the windpop variable for each grid point using do loops..
In order to plot the windpop values, I need to convert the windpop from 1-d
to 3-d (nlev,nlat,mlon)..
I use the onedtond function and the new variable wp has the right dimensions
(1x89x184). However, the values did not correspond with the right values..
The first value which belongs to first grid point (0,0,0) was repeated for
16376 times (up to (0,89,184) grid point), then the same is applied to the
second value etc.
The desired results would be (0,0,0) 1st value, (0,1,0) 2nd etc
I can't understand where is the error..
Has anybody got any idea?
Any help would be appreciated
The source files are available at
<http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1961-19
70_wss.nc.gz>
http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1961-197
0_wss.nc.gz
<http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1971-19
80_wss.nc.gz>
http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1971-198
0_wss.nc.gz
<http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1981-19
90_wss.nc.gz>
http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1981-199
0_wss.nc.gz
for free download.
Attached script
;*************************************************
; Plot wind power potential over Mediterranean
;*************************************************
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"
begin
;=======================================================
; Variables definition and files handling
;=======================================================
DataDir = "/data3/koletsis/Current/A2/"
Files = systemfunc ( " ls " + DataDir + \
"C4IRCA3_A2_ECHAM5_DM_*.nc")
numFiles = dimsizes(Files)
print(Files)
print(" ")
f1 = addfiles(Files,"r")
time = f1[:]->time
ws = f1[:]->wss
ws@_FillValue = 1.e+30
lat2d = f1[0]->lat
lon2d = f1[0]->lon
pollat = 39.25
pollon = -162.0
;======================================================
; Choose time period
;======================================================
time2 = ut_calendar(time, 0)
time3 = ut_calendar(time, -2) ; Date format
year = tointeger(time2(:,0)) ; Convert to integer
month = tointeger(time2(:,1))
day = tointeger(time2(:,2))
beg_year = 1961
end_year = 1990
beg_month = 01
end_month = 12
beg_day = 01
end_day = 31
ymdStrt = beg_year*10000+beg_month*100+beg_day
ymdLast = end_year*10000+end_month*100+end_day
print("ymdStrt="+ymdStrt+" ymdLast="+ymdLast)
itime = ind(time3.ge.ymdStrt.and.time3.le.ymdLast)
if (ismissing(itime(0))) then ; Error message
print("itime is missing:something is going wrong")
exit
end if
;======================================================================
; Select the subregion (Mediterranean Area)
;======================================================================
latS = 30
latN = 46
lonW = -10
lonE = 41
ji = region_ind(lat2d,lon2d,latS,latN,lonW,lonE)
jStrt = ji(0)
jLast = ji(1)
iStrt = ji(2)
iLast = ji(3)
llat2d = lat2d(jStrt:jLast,iStrt:iLast)
llon2d = lon2d(jStrt:jLast,iStrt:iLast)
ws2 = ws(itime,:,jStrt:jLast,iStrt:iLast)
ws1 = 1.2697*ws2 ; wind speed at hub height
dimx = dimsizes(ws1)
nlev = dimx(1)
nlat = dimx(2)
mlon = dimx(3)
; npts = nlat*mlon
; wpower = new((/nlev,nlat,mlon/),"string")
npt = -1
do kl = 0,nlev-1
do nl = 0,nlat-1
do ml = 0,mlon-1
data = ws1(:,kl,nl,ml)
;======================================================================
; Calculate probability function
;======================================================================
opt = True
opt@bin_spacing = 0.2
max_ws1 = max(data)
min_ws1 = 0.
nbins = floattointeger(((max_ws1 - min_ws1) / 0.2)+1) ; bins number
wspdf = pdfx(data,nbins,opt)
; printVarSummary(wspdf)
;=======================================================================
; Calculate cumulative density function
;=======================================================================
xcumsum = cumsum(wspdf,2)
xcumsum_Div = 99.99999 ; in order to avoid
xcumsum = where(xcumsum.ne.100,xcumsum,xcumsum_Div) ; impossible division
;=======================================================================
; Calculate Weibull parameters using least square method
;=======================================================================
yi = log(log(1.0/(1.0 - (xcumsum/100))))
if(any(isnan_ieee(yi))) then ; in order to avoid nan result
value = 2.230200812934578 ; in yi estimation
replace_ieeenan(yi,value,0)
yi@_FillValue = value
end if
yi_Zero = 0
yi = where(yi.ne."-inf",yi,yi_Zero) ; set missing (zero) value -inf
umaxiall = (/0.2,0.4,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,\
3.,3.2,3.4,3.6,3.8,4.,4.2,4.4,4.6,4.8,5.,5.2,5.4,5.6,5.8,\
6.,6.2,6.4,6.6,6.8,7.,7.2,7.4,7.6,7.8,8.,8.2,8.4,8.6,8.8,\
9.,9.2,9.4,9.6,9.8,10.,10.2,10.4,10.6,10.8,11.,11.2,11.4,\
11.6,11.8,12.,12.2,12.4,12.6,12.8,13.,13.2,13.4,13.6,13.8,\
14.,14.2,14.4,14.6,14.8,15.,15.2,15.4,15.6,15.8,16.,16.2,\
16.4,16.6,16.8,17.,17.2,17.4,17.6,17.8,18.,18.2,18.4,18.6,\
18.8,19.,19.2,19.4,19.6,19.8,20.,20.2,20.4,20.6,20.8,21.,\
21.2,21.4,21.6,21.8,22.,22.2,22.4,22.6,22.8,23.,23.2,23.4,\
23.6,23.8,24.,24.2,24.4,24.6,24.8,25.,25.2,25.4,25.6,25.8,\
26.,26.2,26.4,26.6,26.8,27.,27.2,27.4,27.6,27.8,28.,28.2,28.4,\
28.6,28.8,29.,29.2,29.4,29.6,29.8,30.,30.2,30.4,30.6,30.8,\
31.,31.2,31.4,31.6,31.8,32.,32.2,32.4,32.6,32.8,33.,33.2,\
33.4,33.6,33.8,34.,34.2,34.4,34.6,34.8,35.,35.2,35.4,35.6,\
35.8,36.,36.2,36.4,36.6,36.8,37.,37.2,37.4,37.6,37.8,38.,\
38.2,38.4,38.6,38.8,39.,39.2,39.4,39.6,39.8,40./)
umaxi = umaxiall(0:nbins-1)
xi = log(umaxi)
xiyi= xi*yi
sumxy = sum(xiyi(0:nbins-1))
sumy = sum(yi(0:nbins-1))
sumx = sum(xi(0:nbins-1))
xi2 = xi*xi
sumx2 = sum(xi2(0:nbins-1))
b1 = sumy * sumx2
b2 = sumx * sumxy
b3 = nbins * sumx2
b4 = sumx * sumx
b = (b1 - b2) / (b3 - b4)
k1 = nbins * sumxy
k2 = sumx * sumy
k3 = nbins * sumx2
k4 = b4
;========================================================================
; Weibull Shape Parameter k
;========================================================================
k = (k1 - k2) / (k3 - k4)
;========================================================================
; Weibull Scale Parameter c
;========================================================================
ratio = - (b / k)
c = exp(ratio)
;========================================================================
; Calculate wind speed probability density function from 5 up to 25 m/s
;========================================================================
df1 = (k / c)
if(nbins.le.125) then
df2 = (umaxiall(25:nbins) / c) ^ (k-1)
df3 = -((umaxiall(25:nbins) / c) ^ k)
else
df2 = (umaxiall(25:125) / c) ^ (k-1)
df3 = -((umaxiall(25:125) / c) ^ k)
end if
df4 = exp(df3)
df = df1 * df2 * df4
numdf = (dimsizes(df)) - 1
;========================================================================
; Calculate Wind Power Potential
;========================================================================
ginol5a = umaxiall(25) ^ 3 ; epeidh epilegw thn akrh toy bin
if(nbins.le.125) then
ginol25a = umaxiall(nbins) ^ 3 ; to oloklirwma ypologizetai basei
else
ginol25a = umaxiall(125) ^ 3
end if
ginol5b = df(0) * 0.2 ; methodoy toy Simpson
ginol25b = df(numdf) * 0.2
ginol5 = ginol5a * ginol5b
ginol25 = ginol25a * ginol25b
if(nbins.le.125) then
ginol1 = (umaxiall(25:nbins))^3
else
ginol1 = (umaxiall(25:125))^3
end if
ginol2 = df * 0.2
ginol = ginol1 * ginol2
sumol0 = sum(ginol)
sumol = sumol0 + (0.5 * ginol25) - ginol5
windpop = 0.5 * 1.225 * sumol
; npt = npt+1
wp = onedtond(windpop,(/nlev,nlat,mlon/))
; printVarSummary(wp)
delete([/wspdf,data,xcumsum,umaxi,yi,xi,xiyi,xi2,df2,df3,df4,\
df,ginol1,ginol2,ginol/])
end do
end do
end do
end
<test_plot_wp2.ncl>_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
<http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_____
No virus found in this message.
Checked by AVG - www.avg.com
Version: 2013.0.2904 / Virus Database: 3209/6522 - Release Date: 07/26/13
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Jul 29 04:21:15 2013
This archive was generated by hypermail 2.1.8 : Thu Aug 01 2013 - 15:55:03 MDT