# Re: PDF calculation for a grid node

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Fri Jul 12 2013 - 16:35:16 MDT

I'm not sure what you mean by "grid node" versus "grid point".

Did you try callilng the function without any loops?

> wspdf=(pdfx(ws1,nbins,opt))

You may want "ws2" to be a 2D array rather than a 3D array, if you are trying to get the prob density across a lat/lon area for each level:

llat2d = lat2d(jStrt:jLast,iStrt:iLast)
llon2d = lon2d(jStrt:jLast,iStrt:iLast)
nlat = dimsizes(lat2d(:,0))
nlon = dimsizes(lat2d(0,:))
nlev = dimsizes(ws(0,:,0,0))
ws2 = reshape(ws(itime,:,jStrt:jLast,iStrt:iLast),(/nlev,nlat*nlon/))
ws1 = 1.2697*ws2

Then I think you can call pdfx without having to loop.

--Mary

On Jul 11, 2013, at 8:11 AM, Ioannis Koletsis <koletsis@noa.gr> wrote:

> Dear ncl users,
>
> I would like to calculate the pdf of wind speed for a grid node…
> When I calculated it for a grid point everything went well.
> The number of bins, and the maximum wind speed for each grid point have been estimated correctly (line number 92, 89 respectively), for each grid point throughout the time period (10957 records)….
> Moreover, the printVarSummary for max_ws1 and nbins give me the dimensions….
> Number of dimensions 3 (1 x 19 x 17) -> (height x rlat x rlon)
> As you can see at the attached script, I tried to solve the problem using do loops but unfortunately without success…
>
> Has anyone got any idea, how to calculate the pdf for each grid point….
> Any help would be appreciated….
>
> Ioannis
>
>
>
> ;*************************************************
> ; PDF Calculation
> ;*************************************************
>
> begin
>
> 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
>
> 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
> latS = 36
> latN = 40
> lonW = 22
> lonE = 26
>
> 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
>
> ; print(ws1)
> ; printVarSummary(ws2)
> ; printVarSummary(ws1)
>
> ;======================================================================
> ; Calculate probability function
> ;======================================================================
>
> opt = True
> opt@bin_spacing = 0.2
>
> max_ws1 = dim_max_n(ws1,0)
> min_ws1 = 0.
>
> nbins = floattointeger(((max_ws1 - min_ws1) / 0.2)+1) ; bins number
>
> ; print(nbins)
>
> dimx = dimsizes(ws1)
> nlev = dimx(0)
> nlat = dimx(1)
> mlon = dimx(2)
>
> do kk = 1, 10957
> do nl = 0, nlat-1
> do ml = 0, mlon-1
> do nt = 0, nlev-1
> wspdf=(pdfx(ws1(kk,nt,nl,ml),nbins,opt))
> end do
> end do
> end do
> end do
> print(wspdf)
>
> end
>
>
>
> <pdf.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 Fri Jul 12 16:35:21 2013

This archive was generated by hypermail 2.1.8 : Fri Jul 12 2013 - 16:37:39 MDT