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

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
Received on Fri May 19 2006 - 12:56:42 MDT