1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
sh sym clock_secs
use "/home/data/socat/SOCAT2_data_table_b900_f2d1_5213.nc"
let siz = `fco2_recomputed,return=isize`
def sym rtol = $1"0.01"
let npts = $2"`siz`"
let parm = $3"2"
def sym tol = ($rtol)
def sym parm = `parm`
set mem/siz=400
! replace the last datum in each trajectory with the bad-value.
! Normalize the Fco2 data
stat fco2_recomputed
let fco2_norm = (fco2_recomputed - ($stat_mean))/($stat_std)
sh sym clock_secs
let/title="`fco2_recomputed,return=title` normalized"/units="`fco2_recomputed,return=units`" \
fco_norm_with_gaps = separate(fco2_norm[i=1:`npts`], rowsize, 0)
let/title="`fco2_recomputed,return=title`"/units="`fco2_recomputed,return=units`" \
fco_with_gaps = separate(fco2_recomputed[i=1:`npts`], rowsize, 0)
let/units="`longitude,return=units`"/title="`longitude,return=title`" \
lon_with_gaps = separate(longitude[i=1:`npts`], rowsize, 1)
let/title="`latitude,return=title`"/units="`latitude,return=units`" lat_with_gaps = \
separate(latitude[i=1:`npts`], rowsize, 0)
let npts = `lat_with_gaps,return=isize`
! Compute and save the distance between successive lon/lat locations on the globe
! Takes about 8 seconds on dunkel.
! sh sym clock_secs
! let/title="distance on sphere"/units="m" ypts_lonlat = dist2(lon_with_gaps, lat_with_gaps, `parm`)
! save/clobber/file="/home/data/ansley/lonlat_dist_on_sphere.nc" ypts_lonlat
! sh sym clock_secs
! can var ypts_lonlat
! pause
use "/home/data/ansley/lonlat_dist_on_sphere.nc"
! Normalize the distance on the globe.
stat ypts_lonlat[d=2]
let dist_norm = (ypts_lonlat[d=2] - (($stat_mean)))/($stat_std)
set data 1
! Sample using the piecewise linear approximation, applied to the
! distance along the globe of the lon-lat path.
let nf = npts
let ntol = 1
let rest = 1
let cont = 1
let decimate_lonlat = piecewise(dist_norm, ($rtol), ntol, rest, cont)
let nout_lonlat = decimate_lonlat[j=1,i=@ngd]
sh sym clock_secs
load decimate_lonlat
sh sym clock_secs
say Sample with Tolerance `($rtol)` on lon/lat distance on sphere: Returned `nout_lonlat` of original `nf` Points, `nout_lonlat/nf`
pause
def sym nout_lonlat = `nout_lonlat`
say `nout_lonlat/nf`
let xsample_on_lonlat = xsequence(decimate_lonlat[j=1,i=1:($nout_lonlat)])
def axis/x=1:($nout_lonlat):1 xsample_axis
let sample_on_lonlat = xsample_on_lonlat[gx=xsample_axis@asn]
let fco2sample_lonlat = samplei(fco_with_gaps[i=1:`npts`], sample_on_lonlat)
let fco2norm_sample_lonlat = samplei(fco_norm_with_gaps[i=1:`npts`], sample_on_lonlat)
let lonsample_lonlat = samplei(lon_with_gaps[i=1:`npts`], sample_on_lonlat)
let latsample_lonlat = samplei(lat_with_gaps[i=1:`npts`], sample_on_lonlat)
! Now apply the piecewise linear approximation to the remaining points, fitting
! lines through the fco2 measurements
def sym rtol = 10*($tol)
let decimate_fco2 = piecewise(fco2norm_sample_lonlat[i=1:($nout_lonlat)], ($rtol), ntol, rest, cont)
sh sym clock_secs
def sym nout_fco2 = `decimate_fco2[j=1,i=@ngd]`
sh sym clock_secs
say Sample again on fco2, with Tolerance `($rtol)`: Returned ($nout_fco2) of previous ($nout_lonlat) Points , `($nout_fco2)/($nout_lonlat)`
pause
let sample_on_both = xsequence(decimate_fco2[j=1,i=1:($nout_lonlat)])
keymark 1
! 1000000 takes a minute to plot, interactive mode
! STYLES:
! /line/fast/miss=blank
! /sym=dot/miss=blank
set win/siz=0.6
g margins 0,0,0,0
go basemap x=-180:180 y=-90:90 20
LET/title="`fco2_recomputed,return=title` sampled"/units="`fco2_recomputed,return=units`" \
fco2sample = samplei(fco2sample_lonlat[i=1:($nout_lonlat)], sample_on_both)
LET/title="`longitude,return=title`"/units="`longitude,return=units`" \
lonsample = samplei(lonsample_lonlat[i=1:($nout_lonlat)], sample_on_both)
LET/title="`latitude,return=title`"/units="`latitude,return=units`" \
latsample = samplei(latsample_lonlat[i=1:($nout_lonlat)], sample_on_both)
def sym plottitle = "`fco2_recomputed,return=title`<nl>($plottitle)"
def sym vlevels = (-inf)(160,260,20)(260,450,10)(440,560,20)(inf)
let ribbon = 1
IF `ribbon EQ 1` THEN
sh sym clock_secs
plot/over/pal=rnb/vs/ribbon/line/miss=blank/nolab/levels=($vlevels)/key lonsample, latsample, fco2sample
sh sym clock_secs
pause
frame/file=sample_ribbon.gif
exit/script
ENDIF ! ribbon
! or
! Plot the twice-sampled data.
! Note it doesn't help things to pre-load lonsample, latsample etc
! when we define vars in GO POLYMARK, the uvars all get purged.
let twice = 1
if `twice EQ 1` THEN
sh sym clock_secs
go polymark poly/over/pal=rnb/nolab/levels=($vlevels)/key \
lonsample latsample fco2sample square 0.15
sh sym clock_secs
! label/nouser `($ppl$xlen)/2`, -0.8, 0, 0, 0.12, \
! @CR2-step decimation. lon/lat distance, parm ($parm): Sample w/ piecwise linear approx,<nl>\
! 1st on path distance, then on fCO2. Start `npts` to ($($nout_fco2)) pts
frame/file=sample_tol($tol)_parm($parm).gif
ELSE
! Plot the data sampled only using path distance.
set mode diag
sh sym clock_secs
go polymark poly/over/pal=rnb/nolab/levels=($vlevels)/key \
lonsample_lonlat latsample_lonlat fco2sample_lonlat square 0.15
sh sym clock_secs
! label/nouser `($ppl$xlen)/2`, -0.8, 0, 0, 0.12, \
! @CR2-step decimation. lon/lat distance, parm ($parm): Sample w/ piecwise linear approx,<nl>\
! only on path distance, (not on fCO2). Start `npts` to ($nout_lonlat) pts
frame/file=sample_tol($tol)_parm($parm)_distance_only.gif
ENDIF
|
