[ncl-talk] Reading and writing data

Appo derbetini appopson4 at gmail.com
Thu Feb 21 03:00:34 MST 2019


Dear Adam,
With the script you send, i am able to read and write data into netcdf
file.
But when i am trying to use cdo to select specific year  with this command

cdo   selyear,2000   NC131.nc     aaa.nc

this error appears
Warning (cdf_set_dim): Inconsistent dimension definition for lat! dimid =
0;  type = 1;  newtype = 2
Warning (cdf_set_dim): Inconsistent dimension definition for lon! dimid =
0;  type = 2;  newtype = 1
cdo selyear (Warning): Year 2000 not found!
cdo selyear (Abort): No timesteps selected!

Here below the script used.

Regards Appo

begin
 yrStrt = 1920
 yrLast = 2014
 file_name   = "Rain_estimate_2.5DegGrid.txt"
 rain    = asciiread(file_name, -1, "string")

lat = fspan(-23.75, 11.25, 15)
lat at units = "degrees_north"
lat at standard_name = "latitude" ;
lat!0 = "lat"
lat&lat = lat
lat at axis = "X"
lat at standard_name = "latitude" ;
lat at long_name = "latitude" ;

;print(lat)

lon = fspan(6.25, 43.75, 16)
lon at units = "degrees_east"
lon at standard_name = "longitude" ;
lon!0 = "lon"
lon&lon = lon
lon at axis = "Y"
lon at standard_name = "longitude" ;
lon at long_name = "longitude" ;



    time = yyyymm_time(yrStrt, yrLast, "integer")

    time at long_name = "time" ;
    time at standard_name = "time" ;
    time at units = "days since 1900-01-01 00:00:0.0"
    time at calendar = "gregorian" ;
    time at axis = "T"



mod_data = new((/(yrLast-yrStrt+1)*12, dimsizes(lat), dimsizes(lon)/),
float)    ;Set up master array that will be filled in.
mod_data!0  = "time"
mod_data&time = time
mod_data!1 = "lat"
mod_data&lat = lat
mod_data!2 = "lon"
mod_data&lon = lon
mod_data at long_name = "precipitation"           ; assign attributes
mod_data at units = "mm" ;
mod_data at missing_value  = mod_data at _FillValue

nrow = dimsizes(rain)


do gg = 0, nrow-2, 100   ; Run a do loop over each 100 row segment
     temp = rain(gg:gg+99)
     coords = str_split(temp(1), " ")
     ;print(coords)     ; examine results of text split
     lat_coord = tofloat(coords(2))
     lon_coord = tofloat(coords(5))
     ii = 0
     do hh = 5,99   ; operate on rows 5-99 of temp
          data = tofloat(str_split(temp(hh), " "))

          mod_data({data(0)*100+1:data(0)*100+12}, {lat_coord},
{lon_coord}) = (/ data(1:) /)

          delete(data)
    end do
     delete(temp)
    ;print("Now done with lat="+lat_coord+", lon="+lon_coord)
end do


;   Save to netcdf
    ;===================================================================
        ntim  = dimsizes(time)                 ; get dimension sizes
    nlat  = dimsizes(lat)
    nlon  = dimsizes(lon)

        diro = "./"                     ; Output directory
        filo = "NC131.nc"             ; Output file
    system("/bin/rm -f " + diro + filo)    ; remove if exists
    fout  = addfile (diro + filo, "c")  ; open output file

    ;===================================================================
    ; explicitly declare file definition mode. Improve efficiency.
    ;===================================================================
        setfileoption(fout,"DefineMode",True)

    ;===================================================================
    ; create global attributes of the file
    ;===================================================================
        fAtt               = True            ; assign file attributes
    fAtt at title         = "NC131"
    fAtt at source_file   =  "Nicholson SE, Klotter D, Dezfuli AK, Zhou L
(2018b) New rainfall datasets for the Congo Basin and surrounding regions.
J Hydrometeorol 19:1379–1396"
    fAtt at Conventions   = "None"
    fAtt at creation_date = systemfunc ("date")
    fileattdef( fout, fAtt )            ; copy file attributes

    ;===================================================================
    ; predefine the coordinate variables and their dimensionality
    ; Note: to get an UNLIMITED record dimension, we set the dimensionality
    ; to -1 (or the actual size) and set the dimension name to True.
    ;===================================================================
        dimNames = (/"time", "lat", "lon"/)
    dimSizes = (/ -1   ,  nlat,  nlon /)
    dimUnlim = (/ True , False, False/)
    filedimdef(fout, dimNames, dimSizes, dimUnlim)

    ;===================================================================
    ; predefine the the dimensionality of the variables to be written out
    ;===================================================================


       filevardef(fout, "time" ,typeof(time),       getvardims(time))
       filevardef(fout, "lat"  ,typeof(lat),
getvardims(lat))
       filevardef(fout, "lon"  ,typeof(lon),
getvardims(lon))
       filevardef(fout, "pr"   ,typeof(mod_data),
getvardims(mod_data))

                                                              ; different
from name on script
    ;===================================================================
    ; Copy attributes associated with each variable to the file
    ; All attributes associated with each variable will be copied.
    ;====================================================================
       filevarattdef(fout,"time" ,time)                    ; copy time
attributes
       filevarattdef(fout,"lat"  ,lat)                     ; copy lat
attributes
       filevarattdef(fout,"lon"  ,lon)                     ; copy lon
attributes
       filevarattdef(fout,"pr"   ,mod_data)                      ; copy PS
attributes

    ;===================================================================
    ; explicitly exit file definition mode. **NOT REQUIRED**
    ;===================================================================
        setfileoption(fout,"DefineMode",False)

    ;===================================================================
    ; output only the data values since the dimensionality and such have
    ; been predefined. The "(/", "/)" syntax tells NCL to only output the
    ; data values to the predefined locations on the file.
    ;====================================================================
       fout->time   = (/time/)
       fout->lat    = (/lat/)
       fout->lon    = (/lon/)
       fout->pr      = (/mod_data/)


end


Le mer. 20 févr. 2019 à 13:42, Appo derbetini <appopson4 at gmail.com> a
écrit :

> Hi Adam,
> After multiple checking i found that error is coming from the fact that
> l"at" and "lon" were inverted here
>
> lat = fspan(6.25,43.75,16)
> lon = fspan(-23.75, 11.25,15)
>
> instead of
>
> lat = fspan( -23.75, 11.25,15)
> lon = fspan(6.25,43.75,16)
>
> Best regards
>
> Appo
>
>
> Le mer. 20 févr. 2019 à 12:44, Appo derbetini <appopson4 at gmail.com> a
> écrit :
>
>> Dear Adam,
>> Thank you for your kind help
>> The script proposed by you is very clear and optimized.
>> But there is an error appearing at this line
>> mod_data({data(0)*100+1:data(0)*100+12},{lat_coord},{lon_coord}) = (/
>> data(1:) /)
>>
>> Telling
>>
>> fatal:NclOneDValGetClosestIndex: finish coordinate index out of range,
>> can't continue
>> fatal:Could not obtain coordinate indexes, unable to perform subscript
>> fatal:["Execute.c":8637]:Execute: Error occurred at or near line 48 in
>> file read_NCwrite_NC131_Adam.ncl
>>
>> Since a hour I am looking how to correct it without success
>>
>> Thanks again
>>
>> Appo
>>
>> Le mar. 19 févr. 2019 à 21:54, Adam Phillips <asphilli at ucar.edu> a
>> écrit :
>>
>>> Hi Appo,
>>> The nice aspect of your ascii file is that it is regular: Each grid
>>> point's data is given over 100 rows, and each is monthly running from
>>> 1920-2014. I would strongly recommend developing a new script that operates
>>> on each 100 row chunk of data, and that uses coordinate subscripting to put
>>> the data in the proper place in your final data array. One negative to this
>>> method is the possibility exists that multiple nearby data points get
>>> assigned to the same grid point, but this can be countered by setting your
>>> master array grid appropriately. An example:
>>>
>>> ; Untested! Review the code below.
>>>  yrStrt = 1920
>>>  yrLast = 2014
>>>  file_name   = "Rain_estimate_2.5DegGrid.txt"
>>>  rain    = asciiread(file_name, -1, "string")
>>>
>>> lat = fspan(6.25,43.75,16)
>>> lat at units = "degrees_north"
>>> lat!0 = "lat"
>>> lat&lat = lat
>>> lon = fspan(-23.75, 11.25,15)
>>> lon at units = "degrees_east"
>>> lon!0 = "lon"
>>> lon&lon = lon
>>> time = yyyymm_time(yrStrt,yrLast,"integer")
>>> mod_data =
>>> new((/(yrLast-yrStrt+1)*12,dimsizes(lat),dimsizes(lon)/),float,-999.)
>>> ;Set up master array that will be filled in.
>>> mod_data!0  = "time"
>>> mod_data&time = time
>>> mod_data!1 = "lat"
>>> mod_data&lat = lat
>>> mod_data!2 = "lon"
>>> mod_data&lon = lon
>>>
>>> nrow = numAsciiRow("file_name")
>>>
>>> do gg = 0,nrow-1,100   ; Run a do loop over each 100 row segment
>>>      temp = rain(gg:gg+99)
>>>      coords = str_split(temp(1)," ")
>>>      print(coords)     ; examine results of text split
>>>      lat_coord = tofloat(coords(2))
>>>      lon_coord = tofloat(coords(5))
>>>      do hh = 5,99   ; operate on rows 5-99 of temp
>>>           data = tofloat(str_split(temp(hh)," "))
>>>
>>> mod_data({data(0)*100+1:data(0)*100+12},{lat_coord},{lon_coord}) = (/
>>> data(1:) /)
>>>           delete(data)
>>>     end do
>>>      delete(temp)
>>>     print("Now done with lat="+lat_coord+", lon="+lon_coord)
>>> end do
>>>
>>> The above should at least get you going. If further changes need to be
>>> made to the above you will need to make them. If you have any further
>>> questions let ncl-talk know.
>>> Adam
>>>
>>>
>>> On Tue, Feb 19, 2019 at 8:27 AM Appo derbetini <appopson4 at gmail.com>
>>> wrote:
>>>
>>>>
>>>> Dear NCL users,
>>>>
>>>> Here a script that i am using to read rainfall data.
>>>> Data are given for each grid point.
>>>> But points are not regularly distributed.
>>>> for some grid point there is no data
>>>>
>>>> I want to construct  regular matrix (lat, lon) where point without data
>>>> will be created and filled by missing value (-999.0 in the script).
>>>>
>>>> My script is giving wrong result !!!
>>>>
>>>> Any help will be appreciated
>>>>
>>>>
>>>> Regards
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> ncl-talk at ucar.edu
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>
>>>
>>> --
>>> Adam Phillips
>>> Associate Scientist,  Climate and Global Dynamics Laboratory, NCAR
>>> www.cgd.ucar.edu/staff/asphilli/   303-497-1726
>>>
>>> <http://www.cgd.ucar.edu/staff/asphilli>
>>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190221/5eeaa366/attachment.html>


More information about the ncl-talk mailing list