# Re: Onedtond problem

From: Ioannis Koletsis <koletsis_at_nyahnyahspammersnyahnyah>
Date: Mon Jul 29 2013 - 04:21:01 MDT

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

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

<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