Re: How to read a 3-d ascii grid data file ?

From: Adam Phillips <asphilli_at_nyahnyahspammersnyahnyah>
Date: Fri, 19 May 2006 15:42:13 -0600

As a final thought on the issue of reading semi-complicated ascii files
into NCL, Dennis has written up an example reading in the original GHCN
gridded precip field ascii file directly into NCL:

NCL can directly read the original file:
     grid_prcp_1900-2006.dat

The generic structure of the ascii file is

month1 year1
-----------------------
Block of values 36*72 |
with multiple columns |
-----------------------
month2 year1
-----------------------
Block of values 36*72 |
with multiple columns |
-----------------------
        :
        :
        :
month12 yearLast
-----------------------
Block of values 36*72 |
with multiple columns |
-----------------------

The fortran program indicates that the array PRCP is dimensioned
     (36,72,12,106) => (nlat,mlon,months,nyrs)
In fortran, the leftmost dimension varies fastest,
hence, latitude is the fastest varying, longitude the
second fastest, months the third fastest and year the slowest.

In NCL, the rightmost dimension varies fastest. Hence, this would be
     (106,12,72,36) => (nyrs,months,mlon,nlat)

Looking at the file, missing values
(_FillValue) are indicated by -32768

Looking at the fortran code, latitude starts at -87.5
                              longitude starts at -177.5
                              data must be scaled by 0.01

The following can be used. The key functioms are:
     built-in: asciiread and onedtond

The following can be used for coordinates
     contributed:
http://www.ncl.ucar.edu/Document/Functions/Contributed/latGlobeFo.shtml
http://www.ncl.ucar.edu/Document/Functions/Contributed/lonGlobeFo.shtml

;______________________________________________________

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
   nmos = 12
   yrStrt = 1900
   yrLast = 2005

   nyrs = yrLast-yrStrt+1
   ntim = nyrs*nmos

   nlat = 36
   mlon = 72
   npts = nlat*mlon ; number of pts in each block

   diri = "./"
   fili = "grid_prcp_1900-2006.dat"
   data = asciiread (diri+fili, -1, "float") ; read all values
   data@_FillValue = -32768
   printVarSummary(data)

   prcp = new ( (/nyrs,nmos,mlon,nlat/), "float", data@_FillValue)

                                  ; number of values in 'header'
   nhd = 2 ; [month year]
   nStrt = 0
   do nyr=0,nyrs-1
     do nmo=0,nmos-1
        mon = data(nStrt)
        year = data(nStrt+1)
        nLast= nStrt+npts+1

        prcp(nyr,nmo,:,:) = onedtond( data(nStrt+2:nLast)*0.01 ,
(/mlon,nlat/))

        nStrt = nStrt+npts+2
     end do
   end do
                                 ; contributed.ncl
   lat = latGlobeFo (nlat, "lat", "latitude", "degrees_north")
   lon = lonGlobeFo (mlon, "lon", "longitude","degrees_east")
   lon = lon - 180

   prcp!0 = "year"
   prcp!1 = "month"
   prcp!2 = "lon"
   prcp!3 = "lat"

   prcp&year = ispan(yrStrt,yrLast,1)
   prcp&month = ispan(1,12,1)
   prcp&lon = lon
   prcp&lat = lat

   prcp_at_long_name = "PRCP"
   prcp_at_units = "????"

   printVarSummary(prcp)
   printMinMax(prcp, True)

end

Adam Phillips wrote:
> Hi Lin,
>
> I have coded a script that reads in the ascii file output by your f90
> program, parses the data correctly, and eventually puts the data into an
> array dimensioned (year,month,lat,lon). Finally, the data is written to
> a netCDF file. The script is attached here. There was nothing wrong with
> your ascii file, and NCL had no issues reading it in. The problem with
> numAsciiRow arised because numAsciiRow uses rather inefficient code, so
> for large ascii files it may take a while to return a result. (But it
> will return a result eventually.) Saji had a great suggestion using "wc
> -l" instead, and I believe that that coding will be used for numAsciiRow
> in the next release of NCL.
>
> Please double check the results of this script. (Check that the script
> parsed the data correctly. You can do this by printing out timeseries of
> random grid points and making sure the output matches the data found in
> the original ascii file.)
>
> Note that NCL could have been used to read in the original GHCN ascii
> file. I spent my efforts trying to read in your fortran 90 created ascii
> file.
> Adam
> --------------------------------------------------------------
> Adam Phillips asphilli_at_ucar.edu
> National Center for Atmospheric Research tel: (303) 497-1726
> ESSL/CGD/CAS fax: (303) 497-1333
> P.O. Box 3000
> Boulder, CO 80307-3000 http://www.cgd.ucar.edu/cas/asphilli
>
>
> ------------------------------------------------------------------------
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>
> begin
> syear = 1900
> eyear = 2005
> nmo = 12
> nlat = 36
> nlon = 72
>
> nrowpts = nlat*nlon
> nyr = eyear-syear+1
> nrow = stringtointeger(systemfunc("'wc' -l " + "Grid_PRCP_1900-2005_Ascii.dat"+" | awk '{print $1}'" ))
> print("Number of rows in ascii file = "+nrow)
> a = asciiread("Grid_PRCP_1900-2005_Ascii.dat",(/nrow,15/),"float")
>
> arr = new((/nyr*nmo,nlat,nlon/),"float",-327.68) ; create array (time,lat,lon)
> arr!0 = "time"
> arr!1 = "lat"
> arr&lat = fspan(-87.5,87.5,nlat) ; assign coordinate variables for lat,lon (time not needed)
> arr!2 = "lon"
> arr&lon = fspan(-177.5,177.5,nlon)
>
> do gg = 0,nrowpts-1 ; parse data in an efficient way: grab all timesteps for each
> tlat = a(gg,1) ; data point. Much faster than grabbing one timestep per point
> tlon = a(gg,2) ; for all timesteps.
> arr(:,{tlat},{tlon}) = (/ ndtooned(a(gg::nrowpts,3:)) /)
> print("Done with grid pt "+tlat+"N, "+tlon+"E")
> end do
>
> finarr = new((/nyr,nmo,nlat,nlon/),"float",-327.68) ; create array (year,month,lat,lon)
> finarr!0 = "year"
> finarr&year = ispan(syear,eyear,1)
> finarr!1 = "month"
> finarr&month = ispan(1,nmo,1)
> finarr!2 = "lat"
> finarr&lat = arr&lat
> finarr!3 = "lon"
> finarr&lon = arr&lon
>
> cntr = 0
> do gg = 0,nyr-1
> finarr(gg,:,:,:) = (/ arr(cntr:cntr+11,:,:) /) ; transfer (time,lat,lon) -> (year,month,lat,lon)
> cntr = cntr+12
> print("Done with year = "+finarr&year(gg))
> end do
> delete(arr)
> printVarSummary(finarr)
>
> q = addfile("Grid_PRCP_1900-2005.nc","c")
> q->PRCP = finarr
> end
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
--------------------------------------------------------------
Adam Phillips			             asphilli_at_ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
ESSL/CGD/CAS                               fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000	  http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri May 19 2006 - 15:42:13 MDT

This archive was generated by hypermail 2.2.0 : Tue May 23 2006 - 14:10:46 MDT