|
1 #!/usr/local/bin/perl |
|
2 #====================================================================== |
|
3 # P D 0 2 G R D |
|
4 # doc: Wed Aug 30 11:51:22 2006 |
|
5 # dlm: Wed Dec 6 09:18:49 2017 |
|
6 # (c) 2006 A.M. Thurnherr |
|
7 # uE-Info: 33 21 NIL 0 0 72 2 2 4 NIL ofnI |
|
8 #====================================================================== |
|
9 |
|
10 # make GMT grd files from RDI PD0 file |
|
11 |
|
12 # HISTORY: |
|
13 # Aug 30, 2006: - created at end of GRAVILUCK cruise |
|
14 # Aug 31, 2006: - BUG: ensembles/bins were numbered from 0 |
|
15 # - added -d)imensional coords |
|
16 # Sep 1, 2006: - fiddled with registration |
|
17 # Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] |
|
18 # Jun 18, 2009: - BUG: xysize had been called xyside |
|
19 # Mar 17, 2016: - adapted to new Getopt library |
|
20 # May 19, 2016: - adapted to velBeamToEarth() |
|
21 # Dec 6, 2017: - renamed from RDI2grd |
|
22 |
|
23 # NOTES: |
|
24 # - regular grids only => no dimensional time axis for data collected |
|
25 # in multi-ensemble burst mode! |
|
26 |
|
27 # TODO: |
|
28 # - implement soundspeed corretion |
|
29 # - add temporal pre-averaging |
|
30 |
|
31 use Getopt::Std; |
|
32 $0 =~ m{(.*/)[^/]+}; |
|
33 require "$1RDI_PD0_IO.pl"; |
|
34 require "$1RDI_Coords.pl"; |
|
35 |
|
36 use NetCDF; |
|
37 |
|
38 sub dumpVar($$$$$) |
|
39 { |
|
40 my($var,$units,$long_name,$fname,$dimnum) = @_; |
|
41 |
|
42 my($ncId) = NetCDF::create("$opt_b$var.grd",NetCDF::CLOBBER); |
|
43 NetCDF::setfill($ncId,NetCDF::NOFILL); # NetCDF library bug |
|
44 |
|
45 my($sid) = NetCDF::dimdef($ncId,'side',2); |
|
46 my($aid) = NetCDF::dimdef($ncId,'xysize',($le-$fe+1)*($lastBin-$firstBin+1)); |
|
47 |
|
48 my($xrid) = NetCDF::vardef($ncId,'x_range',NetCDF::DOUBLE,[$sid]); |
|
49 my($yrid) = NetCDF::vardef($ncId,'y_range',NetCDF::DOUBLE,[$sid]); |
|
50 my($zrid) = NetCDF::vardef($ncId,'z_range',NetCDF::DOUBLE,[$sid]); |
|
51 my($spid) = NetCDF::vardef($ncId,'spacing',NetCDF::DOUBLE,[$sid]); |
|
52 my($dmid) = NetCDF::vardef($ncId,'dimension',NetCDF::LONG,[$sid]); |
|
53 |
|
54 my($zvid) = NetCDF::vardef($ncId,'z',NetCDF::FLOAT,[$aid]); |
|
55 |
|
56 NetCDF::attput($ncId,NetCDF::GLOBAL,'title',NetCDF::CHAR,$ARGV[0]); |
|
57 NetCDF::attput($ncId,NetCDF::GLOBAL,'source',NetCDF::CHAR,$usage); |
|
58 |
|
59 NetCDF::attput($ncId,$xrid,'units',NetCDF::CHAR, |
|
60 $opt_d ? 'day number' : 'ensemble number'); |
|
61 NetCDF::attput($ncId,$yrid,'units',NetCDF::CHAR, |
|
62 $opt_d ? 'm' : 'bin number'); |
|
63 NetCDF::attput($ncId,$zrid,'units',NetCDF::CHAR,$units); |
|
64 |
|
65 NetCDF::attput($ncId,$zvid,'long_name',NetCDF::CHAR,$long_name); |
|
66 NetCDF::attput($ncId,$zvid,'scale_factor',NetCDF::DOUBLE,1); |
|
67 NetCDF::attput($ncId,$zvid,'add_offset',NetCDF::DOUBLE,0); |
|
68 NetCDF::attput($ncId,$zvid,'node_offset',NetCDF::LONG,0); |
|
69 |
|
70 NetCDF::endef($ncId); |
|
71 |
|
72 if ($opt_d) { # dimensional grid |
|
73 my($ft) = $dta{ENSEMBLE}[$fe]->{DAYNO}; # x_range(side) |
|
74 my($lt) = $dta{ENSEMBLE}[$le]->{DAYNO}; |
|
75 NetCDF::varput1($ncId,$xrid,0,$ft); |
|
76 NetCDF::varput1($ncId,$xrid,1,$lt); |
|
77 |
|
78 NetCDF::varput1($ncId,$yrid,0, # y_range(side) |
|
79 $dta{DISTANCE_TO_BIN1_CENTER} + $firstBin*$dta{BIN_LENGTH}); |
|
80 NetCDF::varput1($ncId,$yrid,1, |
|
81 $dta{DISTANCE_TO_BIN1_CENTER} + $lastBin*$dta{BIN_LENGTH}); |
|
82 |
|
83 NetCDF::varput1($ncId,$spid,0,($lt-$ft)/($le-$fe)); # spacing(side) |
|
84 NetCDF::varput1($ncId,$spid,1,$dta{BIN_LENGTH}); |
|
85 } else { # non-dim grid |
|
86 NetCDF::varput1($ncId,$xrid,0,$fe+1); # x_range(side) |
|
87 NetCDF::varput1($ncId,$xrid,1,$le+1); |
|
88 NetCDF::varput1($ncId,$yrid,0,$firstBin+1); # y_range(side) |
|
89 NetCDF::varput1($ncId,$yrid,1,$lastBin+1); |
|
90 NetCDF::varput1($ncId,$spid,0,1); # spacing(side) |
|
91 NetCDF::varput1($ncId,$spid,1,1); |
|
92 } |
|
93 |
|
94 NetCDF::varput1($ncId,$dmid,0,$le-$fe+1); # dimension(side) |
|
95 NetCDF::varput1($ncId,$dmid,1,$lastBin-$firstBin+1); |
|
96 |
|
97 my($min) = 9e99; # z(xyside) |
|
98 my($max) = -9e99; |
|
99 my(@data); |
|
100 for (my($b)=$lastBin; $b>=$firstBin; $b--) { |
|
101 for (my($e)=$fe; $e<=$le; $e++) { |
|
102 my($v) = $dta{ENSEMBLE}[$e]->{$fname}[$b][$dimnum]; |
|
103 $v = nan unless defined($v); |
|
104 $min = $v if ($v < $min); |
|
105 $max = $v if ($v > $max); |
|
106 push(@data,$v); |
|
107 } |
|
108 } |
|
109 |
|
110 my(@start) = (0); |
|
111 my(@count) = (scalar(@data)); |
|
112 |
|
113 NetCDF::varput($ncId,$zvid,\@start,\@count,\@data); |
|
114 |
|
115 NetCDF::varput1($ncId,$zrid,0,$min); # z_range(side) |
|
116 NetCDF::varput1($ncId,$zrid,1,$max); |
|
117 |
|
118 NetCDF::close($ncId); |
|
119 } |
|
120 |
|
121 #------ |
|
122 # Usage |
|
123 #------ |
|
124 |
|
125 $usage = "$0 @ARGV"; |
|
126 die("Usage: $0 [-M)agnetic <declination>] [-r)ange <first_ens,last_ens>] " . |
|
127 "[output -b)ase <name>] [-d)imensional coordinates] " . |
|
128 "<PD0 file>\n") |
|
129 unless (&getopts("b:dM:r:") && @ARGV == 1); |
|
130 |
|
131 print(STDERR "WARNING: magnetic declination not set!\n") |
|
132 unless defined($opt_M); |
|
133 |
|
134 unless (defined($opt_b)) { |
|
135 $opt_b = "$ARGV[0]_"; |
|
136 $opt_b =~ m{(.*)\.\d\d\d}; |
|
137 } |
|
138 |
|
139 ($first_ens,$last_ens) = split(',',$opt_r) |
|
140 if defined($opt_r); |
|
141 |
|
142 #---------- |
|
143 # Read Data |
|
144 #---------- |
|
145 |
|
146 print(STDERR "Reading $ARGV[0]..."); |
|
147 readData($ARGV[0],\%dta); |
|
148 printf(STDERR "%d complete ensembles\n",scalar(@{$dta{ENSEMBLE}})); |
|
149 |
|
150 #-------------------------------------------------- |
|
151 # Find Good Ensemble Range & Make Earth Coordinates |
|
152 #-------------------------------------------------- |
|
153 |
|
154 print(STDERR "Checking/transforming data..."); |
|
155 $dta{HEADING_BIAS} = -$opt_M; # magnetic declination |
|
156 |
|
157 if ($dta{BEAM_COORDINATES}) { # coords used |
|
158 $beamCoords = 1; |
|
159 } elsif (!$dta{EARTH_COORDINATES}) { |
|
160 die("$ARGV[0]: only beam and earth coordinates implemented so far\n"); |
|
161 } |
|
162 |
|
163 $lastGoodBin = 0; |
|
164 for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities |
|
165 next if (defined($first_ens) && |
|
166 $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens); |
|
167 $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e |
|
168 unless defined($P{first_ens}); |
|
169 last if (defined($last_ens) && |
|
170 $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens); |
|
171 $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER}; |
|
172 $le = $e; |
|
173 |
|
174 die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") |
|
175 if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4); |
|
176 die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n") |
|
177 if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR}); |
|
178 die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") |
|
179 if ($dta{ENSEMBLE}[$e]->{LOW_GAIN}); |
|
180 |
|
181 for (my($b)=0; $b<$dta{N_BINS}; $b++) { |
|
182 next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]) && |
|
183 defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]) && |
|
184 defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]) && |
|
185 defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])); |
|
186 @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = |
|
187 $beamCoords ? velBeamToEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
|
188 : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}); |
|
189 $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1; |
|
190 $good_vels[$b]++; |
|
191 $lastGoodBin = $b if ($b > $lastGoodBin); |
|
192 $firstGoodEns = $e unless defined($firstGoodEns); |
|
193 $lastGoodEns = $e; |
|
194 } |
|
195 } |
|
196 |
|
197 unless (defined($opt_r)) { |
|
198 $fe = $firstGoodEns; |
|
199 $le = $lastGoodEns; |
|
200 } |
|
201 |
|
202 $firstBin = 0; |
|
203 $lastBin = $lastGoodBin; |
|
204 |
|
205 print(STDERR "\n"); |
|
206 print(STDERR "Start: $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n"); |
|
207 print(STDERR "End : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n"); |
|
208 printf(STDERR "Bins : %d-%d\n",$firstBin+1,$lastBin+1); |
|
209 |
|
210 #----------- |
|
211 # Write Data |
|
212 #----------- |
|
213 |
|
214 print(STDERR "Writing NetCDF files"); |
|
215 &dumpVar('u','[m/s]','eastward velocity' ,'VELOCITY',0); print(STDERR '.'); |
|
216 &dumpVar('v','[m/s]','northward velocity','VELOCITY',1); print(STDERR '.'); |
|
217 &dumpVar('w','[m/s]','vertical velocity' ,'VELOCITY',2); print(STDERR '.'); |
|
218 &dumpVar('e','[m/s]','error velocity' ,'VELOCITY',3); print(STDERR '.'); |
|
219 |
|
220 &dumpVar('ea1','[count]','beam-1 echo amplitude','ECHO_AMPLITUDE',0); print(STDERR '.'); |
|
221 &dumpVar('ea2','[count]','beam-2 echo amplitude','ECHO_AMPLITUDE',1); print(STDERR '.'); |
|
222 &dumpVar('ea3','[count]','beam-3 echo amplitude','ECHO_AMPLITUDE',2); print(STDERR '.'); |
|
223 &dumpVar('ea4','[count]','beam-4 echo amplitude','ECHO_AMPLITUDE',3); print(STDERR '.'); |
|
224 |
|
225 &dumpVar('corr1','[count]','beam-1 correlation','CORRELATION',0); print(STDERR '.'); |
|
226 &dumpVar('corr2','[count]','beam-2 correlation','CORRELATION',1); print(STDERR '.'); |
|
227 &dumpVar('corr3','[count]','beam-3 correlation','CORRELATION',2); print(STDERR '.'); |
|
228 &dumpVar('corr4','[count]','beam-4 correlation','CORRELATION',3); print(STDERR '.'); |
|
229 |
|
230 &dumpVar('pcg1','[count]','beam-1 %-good','PERCENT_GOOD',0); print(STDERR '.'); |
|
231 &dumpVar('pcg2','[count]','beam-2 %-good','PERCENT_GOOD',1); print(STDERR '.'); |
|
232 &dumpVar('pcg3','[count]','beam-3 %-good','PERCENT_GOOD',2); print(STDERR '.'); |
|
233 &dumpVar('pcg4','[count]','beam-4 %-good','PERCENT_GOOD',3); print(STDERR '.'); |
|
234 print(STDERR "\n"); |
|
235 |
|
236 exit(0); |