[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