New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
read_sla.pro in branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/dataplot – NEMO

source: branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/dataplot/read_sla.pro @ 3026

Last change on this file since 3026 was 3026, checked in by djlea, 13 years ago

Remove unneeded sec_to_dt from dataplot. Also add obstools build instructions to documentation.
File size: 8.9 KB
Line 
1;+
2PRO  read_sla, Files, NumObs=NumObs, $
3               Latitudes=Latitudes, Longitudes=Longitudes, $
4               ObsSLA=ObsSLA, ModSLA=ModSLA, SLAQC=SLAQC, $
5               MDT=MDT, Satellites=Satellites, Types=Types, Dates=Dates, rmdi=rmdi, $
6          iobs=iobs, jobs=jobs, mstp=mstp, quiet=quiet, track=track
7;+-----------------------------------------------------------
8;IDL program to read in netcdf files of altimeter data.
9;
10;Author:  D. J. Lea       - Feb 2008
11;
12;------------------------------------------------------------
13rmdi = -99999.
14NumFiles=n_elements(Files)
15;RefDate = '1950-01-01'
16;!DATE_SEPARATOR='-'
17RefDate=JULDAY(1,1,1950,0,0)     ; should be at 0 UTC
18
19ifile2=0
20for ifile = 0, NumFiles-1 do begin
21;------------------------------------------------------------
22;1. Open the file containing the data
23  ncid = ncdf_open(Files(ifile), /nowrite)
24  if ncid lt 0 and not keyword_set(quiet) then print, 'Error opening file '+Files(ifile)
25  if ncid ge 0 then begin
26  if not keyword_set(quiet) then print, 'Opened file '+Files(ifile)+' successfully'
27
28;------------------------------------------------------------
29;2. Read in the dimensions in the file
30  ncinfo = ncdf_inquire(ncid)
31;  var = ncdf_varinq( ncid, "SLA" )
32  if not keyword_set(quiet) then print,'ncinfo.ndims ',ncinfo.ndims
33  ncdf_diminq, ncid, 0, name, NData
34  NData2=0
35  cycles=1
36  if (ncinfo.ndims gt 2) then ncdf_diminq, ncid, 2, name, NData2
37  if (ncinfo.ndims gt 1) then ncdf_diminq, ncid, 1, name, cycles
38
39
40  if not keyword_set(quiet) then print,NData, NData2
41
42 NData=max([NData*cycles, NData2*cycles])
43
44; print,NData, NData2
45
46  if (NData gt 0) then begin
47;------------------------------------------------------------
48;3. Allocate the data arrays and read in the data
49  lons   = dblarr(NData)
50  lats   = dblarr(NData)
51  obsval = fltarr(NData)
52  modval = fltarr(NData)
53  mdtval = fltarr(NData) 
54  bytQC  = bytarr(NData)
55  QC     = fltarr(NData)
56  Dats   = dblarr(NData)
57  trackval  = intarr(NData)
58  Sats   = intarr(NData)
59  strSats= strarr(NData)
60  StrVal = strarr(8)
61;  dts    = replicate(!dt_base, NData)
62  dts = dblarr(NData)
63 
64;  info,lons
65  varid = ncdf_varid(ncid, 'LONGITUDE')
66  scale_factor=1.
67  if (varid eq -1) then begin
68   varid = ncdf_varid(ncid, 'Longitudes')
69        ncdf_attget, ncid, varid, 'scale_factor', scale_factor
70  endif
71  ncdf_varget, ncid, varid, lons
72  lons=lons*scale_factor
73;  info, lons
74  if (cycles gt 1) then begin
75    lonsn=lons
76    for i=1,cycles-1 do begin
77       lonsn=[lonsn,lons]
78    endfor
79    lons=lonsn
80  endif
81
82;  info, lons
83
84;  info,lats     
85  varid = ncdf_varid(ncid, 'LATITUDE')
86  scale_factor=1.
87  if (varid eq -1) then begin
88   varid = ncdf_varid(ncid, 'Latitudes')
89        ncdf_attget, ncid, varid, 'scale_factor', scale_factor
90  endif
91  ncdf_varget, ncid, varid, lats
92  lats=lats*scale_factor
93;  info, lats
94  if (cycles gt 1) then begin
95    latsn=lats
96    for i=1,cycles-1 do begin
97       latsn=[latsn,lats]
98    endfor
99    lats=latsn
100  endif
101 
102;  info, lats
103   
104  varid = ncdf_varid(ncid, 'JULD')
105  if (varid ne -1) then begin
106   ncdf_varget, ncid, varid, Dats
107  endif else begin
108        varid = ncdf_varid(ncid, 'BeginDates')
109        ncdf_varget, ncid, varid, Bds
110;        print,bds
111;        info, bds
112        if (cycles gt 1) then begin
113           bds=transpose(bds)
114        endif
115        varid = ncdf_varid(ncid, 'NbPoints')
116        ncdf_varget, ncid, varid, nbpoints
117
118; Read in dataindexes and deltat in cls format file
119; These are used to adjust dates
120        varid = ncdf_varid(ncid, 'DeltaT')
121        if (varid ne -1) then begin
122     ncdf_varget, ncid, varid, deltat
123          ncdf_attget, ncid, varid, 'scale_factor', scale_factor
124          deltat=deltat*scale_factor
125          varid = ncdf_varid(ncid, 'DataIndexes')
126          ncdf_varget, ncid, varid, dataindexes
127         
128        endif
129
130   print, 'deltat ',deltat
131;        print, 'dataindexes ',dataindexes
132
133;        info, nbpoints
134;        print,nbpoints
135;        info, cycles
136        cumbds=[0]
137        nbpointst=nbpoints
138        dataindexest=dataindexes
139        for i=1,cycles-1 do begin
140      nbpointst=[nbpointst,nbpoints]
141                dataindexest=[dataindexest,dataindexes]
142        endfor
143        cumbds=[0,total(nbpointst,/cumulative)]
144;       if not keyword_set(quiet) then print,'cumbds ',cumbds
145        nel=n_elements(Bds)
146;        info,nel
147;       if not keyword_set(quiet) then print,n_elements(dats)
148        for i=0,nel-1 do begin
149;        print,'*',i, cumbds(i+1)-cumbds(i)
150;           print,i,cumbds(i),cumbds(i+1)-1
151           dats(cumbds(i):cumbds(i+1)-1)=bds(i)
152           trackval(cumbds(i):cumbds(i+1)-1)=i
153           dataindexest(cumbds(i):cumbds(i+1)-1)=dataindexest(cumbds(i):cumbds(i+1)-1)-dataindexest(cumbds(i))
154        endfor
155       
156; adjust dats based on dataindex
157   dats=dats+dataindexest*deltat/86400.
158       
159       
160;        print, dats
161  endelse
162
163
164         
165  dts=RefDate + Dats
166   
167  varid = ncdf_varid(ncid, 'MISSION')
168  if (varid ne -1) then begin
169  ncdf_varget, ncid, varid, Sats
170  ncdf_attget, ncid, varid, 'Value_0', StrVal0           
171  ncdf_attget, ncid, varid, 'Value_1', StrVal1 
172  ncdf_attget, ncid, varid, 'Value_2', StrVal2 
173  ncdf_attget, ncid, varid, 'Value_3', StrVal3     
174  ncdf_attget, ncid, varid, 'Value_4', StrVal4 
175  ncdf_attget, ncid, varid, 'Value_5', StrVal5
176  ncdf_attget, ncid, varid, 'Value_6', StrVal6   
177  ncdf_attget, ncid, varid, 'Value_7', StrVal7
178  StrVal=[StrVal0, StrVal1, StrVal2, StrVal3, StrVal4, $
179   StrVal5, StrVal6, StrVal7]   
180  for ival = 0, 7 do begin
181    pts = where(Sats eq ival, count)
182    if count gt 0 then strSats(pts) = StrVal(ival)
183  endfor
184  endif
185
186  varid = ncdf_varid(ncid, 'SLA')
187  ncdf_varget, ncid, varid, obsval
188  ncdf_attget, ncid, varid, '_FillValue', FillValue
189  scale_factor=1.
190  if (ncinfo.ndims gt 2) then ncdf_attget, ncid, varid, 'scale_factor', scale_factor
191  if not keyword_set(quiet) then print,'SLA scale factor: ',scale_factor
192  if not keyword_set(quiet) then print,'SLA FillValue: ',FillValue
193  pts = where(obsval eq FillValue, count)
194;  if not keyword_set(quiet) then print,count
195  obsval=obsval*scale_factor
196  if count gt 0 then obsval(pts) = rmdi
197 
198  varid = ncdf_varid(ncid, 'SLA_Hx')
199  if (varid ne -1) then begin
200  ncdf_varget, ncid, varid, modval
201  ncdf_attget, ncid, varid, '_FillValue', FillValue
202;  if not keyword_set(quiet) then print,'SLA_Hx FillValue ',FillValue
203  pts = where(modval eq FillValue or modval eq -9999.0, count)
204;  if not keyword_set(quiet) then print,count
205  if count gt 0 then modval(pts) = rmdi
206  endif
207   
208  varid = ncdf_varid(ncid, 'SLA_QC')
209  if (varid ne -1) then begin
210  ncdf_varget, ncid, varid, bytQC
211  ncdf_attget, ncid, varid, '_FillValue', FillValue
212  QC(*) = 0.
213  pts = where(bytQC eq FillValue, count)
214  if count gt 0 then QC(pts) = rmdi
215  pts = where(bytQC ne 48, count)
216  if count gt 0 then QC(pts) = 1.
217  endif
218 
219  varid = ncdf_varid(ncid, 'MDT')
220  if (varid ne -1) then ncdf_varget, ncid, varid, mdtval
221
222  varid = ncdf_varid(ncid, 'data_source')
223  if varid ne -1 then begin
224     ncdf_varget, ncid, varid, obstyp
225  endif else begin
226     varid = ncdf_varid(ncid, 'MISSION')
227     if varid ne -1 then ncdf_varget, ncid, varid, obstyp
228  endelse
229
230; reads in multi cycle data the wrong way round!
231  if (cycles gt 1) then begin
232      obsval=reform(obsval,[cycles,ndata/cycles])
233      obsval=transpose(obsval)
234      modval=reform(modval,[cycles,ndata/cycles])
235      modval=transpose(modval)
236      QC=reform(QC,[cycles,ndata/cycles])
237      QC=transpose(QC)
238  endif
239
240  varid = ncdf_varid(ncid, 'IOBS')
241  if (varid ne -1) then ncdf_varget, ncid, varid, iobsval
242
243  varid = ncdf_varid(ncid, 'JOBS')
244  if (varid ne -1) then ncdf_varget, ncid, varid, jobsval
245
246  varid = ncdf_varid(ncid, 'MSTP')
247  if (varid ne -1) then ncdf_varget, ncid, varid, mstpval
248
249;  info, obssla
250;  info, obsval
251;  info, modsla
252;  info, modval
253 
254  if ifile2 eq 0 then begin
255    Latitudes = [float(lats)]
256    Longitudes = [float(lons)]
257    ObsSLA = [obsval(*)]
258    ModSLA = [modval(*)]
259    MDT = [mdtval]
260    SLAQC = [QC]
261    Dates = [dts]
262    Satellites = [strSats]
263    if (n_elements(obstyp) gt 0) then Types = [obstyp]
264    if (n_elements(iobsval) gt 0) then iobs = [iobsval]
265    if (n_elements(jobsval) gt 0) then jobs = [jobsval]
266   if (n_elements(mstpval) gt 0) then mstp = [mstpval]
267   if (n_elements(trackval) gt 0) then track = [trackval]
268  endif else begin
269    Latitudes = [Latitudes, float(lats)]
270    Longitudes = [Longitudes, float(lons)]
271    ObsSLA = [ObsSLA, obsval(*)]
272    ModSLA = [ModSLA, modval(*)]
273    MDT = [MDT, mdtval]   
274    SLAQC = [SLAQC, QC]
275    Dates = [Dates, dts]   
276    Satellites = [Satellites, strSats]
277    if (n_elements(obstyp) gt 0) then Types = [Types, obstyp]
278    if (n_elements(iobsval) gt 0) then iobs = [iobs, iobsval]
279    if (n_elements(jobsval) gt 0) then jobs = [jobs, jobsval]
280    if (n_elements(mstpval) gt 0) then mstp = [mstp, mstpval]
281    if (n_elements(trackval) gt 0) then tracks = [track, trackval]
282  endelse     
283
284ifile2=ifile2+1
285endif
286  ncdf_close, ncid
287
288endif
289
290endfor ; ifile 
291 
292NumObs = n_elements(Latitudes)
293
294END
Note: See TracBrowser for help on using the repository browser.