Re: PDF calculation for a grid node

Dear Mary,

first of all, thanks for your reply..

Unfortunately, I cannot use function reshape because of my ncl version
(6.0.0-beta) as the reshape can be used only in version 6.1.0 and later.

However, I try to reshape my variable using statement new, according to the
attached script..

But, again the calculation of pdf for each grid point of my area has been
failed..

The following message was appeared:

Number of elements of dimension (0) of argument (1) is (323) in function
(pdfx), expected (1) elements..

It seems that pdfx function works only for 1-D variable, but however I
wonder if there is any way to calculate the pdf for each grid point in a
selected area?

Have you got any other idea?

Thanks in advance

Ioannis

Attached Script

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

; Calculation of PDF

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

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

  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

  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)

  nlat = dimsizes(llat2d(:,0))

  nlon = dimsizes(llon2d(0,:))

  nlev = dimsizes(ws(0,:,0,0))

  ws3 = ws(itime,:,jStrt:jLast,iStrt:iLast)

  ws2 = new((/nlev,nlat*nlon/),typeof(ws3))

  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

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

  print(wspdf)

  end

From: ncl-talk-bounces@ucar.edu [mailto:ncl-talk-bounces@ucar.edu] On Behalf
Of Mary Haley
Sent: Saturday, July 13, 2013 1:35 AM
To: Ioannis Koletsis
Cc: ncl-talk@ucar.edu
Subject: Re: PDF calculation for a grid node

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

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

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

  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

  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>
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: 3204/6490 - Release Date: 07/14/13

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