[ncl-talk] Reading and writing data
Dennis Shea
shea at ucar.edu
Thu Feb 21 08:20:38 MST 2019
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/976c2cb5/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: read_NCwrite_NC131.ncl
Type: application/octet-stream
Size: 6564 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190221/976c2cb5/attachment.obj>
More information about the ncl-talk
mailing list