Re: Onedtond problem

From: Rick Brownrigg <brownrig_at_nyahnyahspammersnyahnyah>
Date: Thu Jul 25 2013 - 09:45:24 MDT

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-1970_wss.nc.gz
> http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1971-1980_wss.nc.gz
> http://ensemblesrt3.dmi.dk/data/A2/C4I/DM/C4IRCA3_A2_ECHAM5_DM_25km_1981-1990_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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jul 25 09:45:43 2013

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