Re: Re: [ncl-talk] How to read a 3-d ascii grid data file ?

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Thu, 18 May 2006 11:01:26 -0600 (MDT)

>Thanks for the reply. I have set my data as the format below:
>
>190001.0 87.50 -157.50 -327.68 -327.68 -327.68 -327.68 -327.68
-327.68 -327.68
>190001.0 87.50 -152.50 -327.68 -327.68 -327.68 -327.68 -327.68
-327.68 -327.68
>........
>the 1st column is yyyymm (float) and 2nd/3rd is lat/lon, the other right 12
colums is the value of Jan to Dec.
>
>I use the code below to read it,But still can't get the right result.
>(NCL use more than 95% of my xeon cpu£¨but fail to finish running£©
>
>
>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
>filName = "Grid_PRCP_1900-2005_Ascii.dat"
>nRow = numAsciiRow (filName)
>nCol = 15
>
>data = asciiread (filName, (/nRow,nCol/), "float")
>printVarSummary(data) ; It seem that ncl stop here because the
printVarSummary(data)
> didn't give a output
>
>yyyymm = floattointeger(data(:,0)) ; this convert float to integer
>
>lats = data(0:2591,1) ; 36*72=2592 grids
>lons = data(0:2591,2)
>PRCP = data(:,3:)
>
>PRCP!0 = 'lat'
>PRCP!1 = 'lon'
>PRCP!2 = 'time'
>PRCP_at_lat = lats
>PRCP_at_lon = lons
>PRCP_at_time = ? ;I also want to know how to set time property of
PRCP appropriately.
________________________________________

Well, first of all, PRCP is not a lat/lon grid.

It is PRCP(nROW,12) where the 12 represents months of the year.
I suggest that you use "printVarSummary" often to look
at an overall view of your variable.

The following is untested and should give you a start.

begin
 filName = "Grid_PRCP_1900-2005_Ascii.dat"
 nRow = numAsciiRow (filName)
 print("nRow="+nRow)
 nCol = 15
 
 data = asciiread (filName, (/nRow,nCol/), "float")
 printVarSummary(data) ; It seem that ncl stop here because the
                                          didn't give a output
 

 DATA = data(:,3:) ; (nRow,12)
 
 yrStrt = 1900
 yrLast = 2005
 nyrs = yrLast-yrStrt+1
 nmos = 12
 
 ntim = nyrs*nmos
 time = new ( ntim, "integer")
 nt = 0
 do yr=yrStrt,yrLast
    time(nt:nt+nmos-1) = yr*100 + ispan(1,nmos,1)
    nt = nt+nmos
 end do
 time_at_long_name = "year_month"
 time_at_units = "yyyymm"
 time!0 = "time"
 time&time = time
 printVarSummary(time)
 
 nlat =
 mlon =

 latStrt = 87.5
 lonStrt = -157.5
 
 dlat = 2.5 ; ???
 dlon = 2.5
 
 lon = ispan(0,mlon-1)*dlon + lonStrt
 lon!0 = "lon"
 lon_at_long_name = "longitude"
 lon_at_units = "degrees_east"
 lon&lon = lon
 printVarSummary(lon)
 
  
 lat = ispan(0,mlon-1)*dlon + lonStrt
 lat!0 = "lon"
 lat_at_long_name = "latitude"
 lat_at_units = "degrees_north"
 lat&lat = lat
 printVarSummary(lat)
 
 PRCP = new ( (/ntim,nlat,mlon/), "float")
 PRCP!0 = "time"
 PRCP!1 = "lat"
 PRCP!2 = "lon"
 PRCP&lat = lat
 PRCP&lon = lon
 PRCP&time = time
 PRCP_at_long_name = "..."
 PRCP_at_units = "..."
 
 nt = 0
 do nr=0,nRow-1 ; loop over rows
    PRCP(nt:nt+11,:,:) = (/ DATA(nr,:) /) ; store 12 months of data
    nt = nt+12
 end do
 
 printVarSummary(PRCP)
 printMinMax(PRCP, True)
 
end
 
Please do not email me directly. I do see ncl-talk

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu May 18 2006 - 11:01:26 MDT

This archive was generated by hypermail 2.2.0 : Thu May 18 2006 - 15:57:32 MDT