[ncl-talk] Reading and writing data
Appo derbetini
appopson4 at gmail.com
Thu Feb 21 08:50:07 MST 2019
Thank you very much
The script is nicely tuned now
regards
Appo
Le jeu. 21 févr. 2019 à 16:20, Dennis Shea <shea at ucar.edu> a écrit :
> CDO is not an NCAR product. Hence, not much experience with its
> expectations.
>
> I speculate that *time* should be something like: "days since
> 1900-01-01 00:00:00"
>
> ;---Define time period
>
> yrStrt = 1920
> yrLast = 2014
> yyyymm = *yyyymm_time*
> <https://www.ncl.ucar.edu/Document/Functions/Contributed/yyyymm_time.shtml>(yrStrt,yrLast,"integer")
> ; yyyymm(time)
> ntim = *dimsizes*(yyyymm) ; all years
> & months
> *printVarSummary*(yyyymm)
> *print*("---")
>
> ; Create time for netCDF
>
> yyyy = yyyymm/100
> mm = yyyymm - yyyy*100
> dd = *conform*
> <https://www.ncl.ucar.edu/Document/Functions/Built-in/conform.shtml>(yyyy,
> 1, -1)
> hh = *conform*(yyyy, 0, -1)
> mn = *conform*(yyyy, 0, -1)
> sc = *conform*(yyyy, 0, -1)
> tunits = "days since 1900-01-01 00:00:00"
>
> time = *cd_inv_calendar*
> <https://www.ncl.ucar.edu/Document/Functions/Built-in/cd_inv_calendar.shtml>(yyyy,mm,dd,hh,mn,sc,tunits,
> 0)
> time!0 = "time"
> time&time = time
> time at units= tunits
>
> yyyymm!0 = "time"
> yyyymm&time := time
>
> See attached
>
> D
>
> On Thu, Feb 21, 2019 at 3:01 AM Appo derbetini <appopson4 at gmail.com>
> wrote:
>
>> 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>
>>>>>
>>>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190221/06cb2fa6/attachment.html>
More information about the ncl-talk
mailing list