# Onedtond problem

From: Ioannis Koletsis <koletsis_at_nyahnyahspammersnyahnyah>
Date: Thu Jul 25 2013 - 07:17:57 MDT

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

Attached script

;*************************************************

; Plot wind power potential over Mediterranean

;*************************************************

begin

;=======================================================

; Variables definition and files handling

;=======================================================

Files = systemfunc ( " ls " + DataDir + \

"C4IRCA3_A2_ECHAM5_DM_*.nc")

numFiles = dimsizes(Files)

print(Files)

print(" ")

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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Thu Jul 25 07:18:22 2013

This archive was generated by hypermail 2.1.8 : Thu Jul 25 2013 - 21:02:42 MDT