|
1
|
|
|
# -*- coding: utf-8 -*- |
|
2
|
|
|
# vim:fileencoding=utf-8 |
|
3
|
|
|
# |
|
4
|
|
|
# Copyright (c) 2017 Stefan Bender |
|
5
|
|
|
# |
|
6
|
|
|
# This file is part of sciapy. |
|
7
|
|
|
# sciapy is free software: you can redistribute it or modify it |
|
8
|
|
|
# under the terms of the GNU General Public License as published by |
|
9
|
|
|
# the Free Software Foundation, version 2. |
|
10
|
|
|
# See accompanying LICENSE file or http://www.gnu.org/licenses/gpl-2.0.html. |
|
11
|
1 |
|
"""SCIAMACHY level 1c limb spectra hdf5 interface |
|
12
|
|
|
""" |
|
13
|
|
|
|
|
14
|
1 |
|
from __future__ import absolute_import, division, print_function |
|
15
|
|
|
|
|
16
|
1 |
|
import logging |
|
17
|
|
|
|
|
18
|
1 |
|
import numpy as np |
|
19
|
1 |
|
from astropy.time import Time |
|
20
|
|
|
|
|
21
|
1 |
|
from ._types import _limb_data_dtype |
|
22
|
|
|
|
|
23
|
1 |
|
logging.basicConfig(level=logging.INFO, |
|
24
|
|
|
format="[%(levelname)-8s] (%(asctime)s) " |
|
25
|
|
|
"%(filename)s:%(lineno)d %(message)s", |
|
26
|
|
|
datefmt="%Y-%m-%d %H:%M:%S %z") |
|
27
|
|
|
|
|
28
|
1 |
|
def _calc_angles(sza_in, los_in, saa_in, z_in, z_out, earth_radius): |
|
29
|
|
|
"""SCIAMACHY limb scan angle setup |
|
30
|
|
|
|
|
31
|
|
|
Calculates the solar zenith, solar azimuth, and line-of-sight angles |
|
32
|
|
|
to the tangent point (tp), to the top of the atmosphere (toa), |
|
33
|
|
|
and to the satellite (sat). All angles are in degrees. |
|
34
|
|
|
|
|
35
|
|
|
Parameters |
|
36
|
|
|
---------- |
|
37
|
|
|
sza_in : ndarray |
|
38
|
|
|
The solar zenith angles (per tangent point). |
|
39
|
|
|
los_in : ndarray |
|
40
|
|
|
The line-of-sight zenith angles (per tangent point). |
|
41
|
|
|
saa_in : ndarray |
|
42
|
|
|
The relative solar azimuth angles (per tangent point). |
|
43
|
|
|
z_in : ndarray |
|
44
|
|
|
The input altitudes, can be tangent heights, top of |
|
45
|
|
|
atmosphere, or satellite altitudes. |
|
46
|
|
|
z_out : ndarray |
|
47
|
|
|
The output altitude, can be tangent heights, top of |
|
48
|
|
|
atmosphere, or satellite altitudes. |
|
49
|
|
|
earth_radius : ndarray |
|
50
|
|
|
The earth radii (per tangent point). |
|
51
|
|
|
|
|
52
|
|
|
Returns |
|
53
|
|
|
------- |
|
54
|
|
|
sza_out : ndarray |
|
55
|
|
|
The calculated solar zenith angles for the output altitudes. |
|
56
|
|
|
los_out : ndarray |
|
57
|
|
|
The calculated line-of-sight angles for the output altitudes. |
|
58
|
|
|
saa_out : ndarray |
|
59
|
|
|
The calculated solar azimuth angles for the output altitudes. |
|
60
|
|
|
""" |
|
61
|
|
|
def sign(x, y): |
|
62
|
|
|
xa = np.full_like(y, np.abs(x)) |
|
63
|
|
|
xa[np.where(y < 0.)] *= -1 |
|
64
|
|
|
return xa |
|
65
|
|
|
|
|
66
|
|
|
cos_psi_in = np.cos(np.radians(sza_in)) |
|
67
|
|
|
mju_in = np.cos(np.radians(los_in)) |
|
68
|
|
|
cos_phi_0 = np.cos(np.radians(saa_in)) |
|
69
|
|
|
sin_phi_0 = np.sin(np.radians(saa_in)) |
|
70
|
|
|
|
|
71
|
|
|
# /* start original calculation in angles.f */ |
|
72
|
|
|
r = earth_radius + z_in |
|
73
|
|
|
sin_1 = np.sqrt(1.0 - mju_in**2) |
|
74
|
|
|
h_0 = r * (sin_1 - 1.0) + z_in |
|
75
|
|
|
|
|
76
|
|
|
z_out[np.where(z_out < h_0)] = h_0[np.where(z_out < h_0)] |
|
77
|
|
|
|
|
78
|
|
|
delta = (np.sqrt((2.0 * earth_radius + z_in + h_0) * (z_in - h_0)) - |
|
79
|
|
|
np.sqrt((2.0 * earth_radius + z_out + h_0) * (z_out - h_0))) |
|
80
|
|
|
mju_out = ((mju_in * r - delta) / |
|
81
|
|
|
np.sqrt((mju_in * r - delta)**2 + (r * sin_1)**2)) |
|
82
|
|
|
sin_out = r * sin_1 / (earth_radius + z_out) |
|
83
|
|
|
sin_psi = np.sqrt(1.0 - cos_psi_in**2) |
|
84
|
|
|
zeta_0 = mju_in * cos_psi_in - sin_1 * sin_psi * cos_phi_0 |
|
85
|
|
|
ksi_0 = mju_in * sin_psi + sin_1 * cos_psi_in * cos_phi_0 |
|
86
|
|
|
|
|
87
|
|
|
cos_psi_out = ((cos_psi_in * r - delta * zeta_0) / |
|
88
|
|
|
np.sqrt((r * cos_psi_in - delta * zeta_0)**2 + |
|
89
|
|
|
(r * sin_psi - delta * ksi_0)**2 + |
|
90
|
|
|
(delta * sin_1 * sin_phi_0)**2)) |
|
91
|
|
|
sin_psi_out = np.sqrt(1.0 - np.clip(cos_psi_out * cos_psi_out, 1.0, np.inf)) |
|
92
|
|
|
#sin_psi_out = np.sqrt(1.0 - np.clip(cos_psi_out*cos_psi_out, -np.inf, 1.0)); |
|
93
|
|
|
|
|
94
|
|
|
eta_0 = ((r - delta * mju_in) * sin_psi * cos_phi_0 - delta * sin_1 * cos_psi_in) / (earth_radius + z_out) |
|
95
|
|
|
|
|
96
|
|
|
eta_0 = eta_0 + 1.0e-39 # /*! numerical stabilization*/ |
|
97
|
|
|
|
|
98
|
|
|
s1 = eta_0 / np.sqrt(eta_0**2 + (sin_psi * sin_phi_0)**2) |
|
99
|
|
|
s1 = s1 - sign(1.0e-13, s1) # /* ! numerical stabilization*/ |
|
100
|
|
|
|
|
101
|
|
|
sd = r * sin_psi * sin_1 * sin_phi_0 / (earth_radius + z_out) / (sin_psi_out * sin_out + 1.0e-78) |
|
102
|
|
|
|
|
103
|
|
|
phi_out = sign(np.arccos(s1), sd) |
|
104
|
|
|
|
|
105
|
|
|
sza_out = np.degrees(np.arccos(cos_psi_out)) |
|
106
|
|
|
los_out = np.degrees(np.arccos(mju_out)) |
|
107
|
|
|
saa_out = np.degrees(phi_out) |
|
108
|
|
|
# /* set back direction of line of sight |
|
109
|
|
|
# if (z_out > z_in) |
|
110
|
|
|
# *saa_out -= PI; */ |
|
111
|
|
|
logging.debug("calculated tangent_height: %s", h_0) |
|
112
|
|
|
return (sza_out, los_out, saa_out) |
|
113
|
|
|
|
|
114
|
1 |
|
def _middle_coord(lat1, lon1, lat2, lon2): |
|
115
|
|
|
sin_lat = 0.5 * (np.sin(np.radians(lat1)) + np.sin(np.radians(lat2))) |
|
116
|
|
|
cos_lat = 0.5 * (np.cos(np.radians(lat1)) + np.cos(np.radians(lat2))) |
|
117
|
|
|
sin_lon = 0.5 * (np.sin(np.radians(lon1)) + np.sin(np.radians(lon2))) |
|
118
|
|
|
cos_lon = 0.5 * (np.cos(np.radians(lon1)) + np.cos(np.radians(lon2))) |
|
119
|
|
|
return (np.degrees(np.arctan2(sin_lat, cos_lat)), |
|
120
|
|
|
np.degrees(np.arctan2(sin_lon, cos_lon))) |
|
121
|
|
|
|
|
122
|
1 |
|
_state_txt = { |
|
123
|
|
|
27: "SCIAMACHY limb eclips", |
|
124
|
|
|
28: "SCIAMACHY limb", |
|
125
|
|
|
29: "SCIAMACHY limb", |
|
126
|
|
|
30: "SCIAMACHY limb", |
|
127
|
|
|
31: "SCIAMACHY limb", |
|
128
|
|
|
32: "SCIAMACHY limb", |
|
129
|
|
|
33: "SCIAMACHY limb", |
|
130
|
|
|
55: "SCIAMACHY limb mesosp", |
|
131
|
|
|
} |
|
132
|
|
|
|
|
133
|
1 |
|
def read_hdf5_limb_state_common_data(self, hf, lstate_id, state_in_orbit, cl_id): |
|
134
|
|
|
"""SCIAMACHY level 1c common data |
|
135
|
|
|
|
|
136
|
|
|
Parameters |
|
137
|
|
|
---------- |
|
138
|
|
|
hf : opened file |
|
139
|
|
|
Pointer to the opened level 1c HDF5 file |
|
140
|
|
|
lstate_id : int |
|
141
|
|
|
The limb state id. |
|
142
|
|
|
state_in_orbit : int |
|
143
|
|
|
The number in this batch of states for the header. |
|
144
|
|
|
cl_id : int |
|
145
|
|
|
The spectral cluster number. |
|
146
|
|
|
|
|
147
|
|
|
Returns |
|
148
|
|
|
------- |
|
149
|
|
|
success : int |
|
150
|
|
|
0 on success, |
|
151
|
|
|
1 if an error occurred, for example if the measurement data |
|
152
|
|
|
set for the requested limb and cluster ids is empty. |
|
153
|
|
|
""" |
|
154
|
|
|
# MDS = measurement data set |
|
155
|
|
|
cl_mds_group_name = "/MDS/limb_{0:02d}/cluster_{1:02d}".format(lstate_id, cl_id + 1) |
|
156
|
|
|
cl_mds_group = hf.get(cl_mds_group_name) |
|
157
|
|
|
if cl_mds_group is None: |
|
158
|
|
|
return 1 |
|
159
|
|
|
# Load meta data |
|
160
|
|
|
self.metadata["calibration"] = hf.attrs["Calibration"].decode() |
|
161
|
|
|
self.metadata["l1b_product"] = hf.get("/MPH")["product_name"][0].decode() |
|
162
|
|
|
self.metadata["orbit"] = hf.get("/MPH")["abs_orbit"][0] |
|
163
|
|
|
self.metadata["state_id"] = hf.get("/ADS/STATES")["state_id"][lstate_id] |
|
164
|
|
|
self.metadata["software_version"] = hf.get("/MPH")["software_version"][0].decode() |
|
165
|
|
|
self.metadata["keyfile_version"] = hf.get("/SPH")["key_data_version"][0].decode() |
|
166
|
|
|
self.metadata["mfactor_version"] = hf.get("/SPH")["m_factor_version"][0].decode() |
|
167
|
|
|
init_ver = hf.get("/SPH")["init_version"][0].decode().strip() |
|
168
|
|
|
self.metadata["init_version"], decont = init_ver.split() |
|
169
|
|
|
self.metadata["decont_flags"] = decont.lstrip("DECONT=") |
|
170
|
|
|
orb_phase = hf.get("/ADS/STATES")["orb_phase"][lstate_id] |
|
171
|
|
|
j_day_0 = 2451544.5 # 2000-01-01 |
|
172
|
|
|
dsr_d, dsr_s, dsr_us = hf.get("/ADS/STATES")["dsr_time"][lstate_id] |
|
173
|
|
|
state_dt = Time(dsr_d + j_day_0 + dsr_s / 86400. + dsr_us / (86400. * 1e6), |
|
174
|
|
|
format="jd").datetime |
|
175
|
|
|
self.metadata["date"] = state_dt.strftime("%d-%b-%Y %H:%M:%S.%f") |
|
176
|
|
|
|
|
177
|
|
|
logging.debug("applied calibrations: %s", self.metadata["calibration"]) |
|
178
|
|
|
logging.debug("product: %s, orbit_nr: %s, state_id: %s, orb_phase: %s", |
|
179
|
|
|
self.metadata["l1b_product"], self.metadata["orbit"], |
|
180
|
|
|
self.metadata["state_id"], orb_phase) |
|
181
|
|
|
logging.debug("soft_ver: %s, key_ver: %s, mf_ver: %s, init_ver: %s, " |
|
182
|
|
|
"decont_ver: %s", self.metadata["software_version"], |
|
183
|
|
|
self.metadata["keyfile_version"], self.metadata["mfactor_version"], |
|
184
|
|
|
self.metadata["init_version"], self.metadata["decont_flags"]) |
|
185
|
|
|
|
|
186
|
|
|
ads_state = hf.get("/ADS/STATES")[lstate_id] |
|
187
|
|
|
cl_n_readouts = ads_state["clus_config"]["num_readouts"][cl_id] |
|
188
|
|
|
cl_intg_time = ads_state["clus_config"]["intg_time"][cl_id] |
|
189
|
|
|
self.metadata["nr_profile"] = 24 // (cl_intg_time * cl_n_readouts) |
|
190
|
|
|
self.metadata["act_profile"] = 0 # always zero for now |
|
191
|
|
|
|
|
192
|
|
|
# Prepare the header |
|
193
|
|
|
try: |
|
194
|
|
|
self.metadata["datatype_txt"] = _state_txt[self.metadata["state_id"]] |
|
195
|
|
|
except KeyError: |
|
196
|
|
|
logging.warn("State id %s not supported.", self.metadata["state_id"]) |
|
197
|
|
|
return 1 |
|
198
|
|
|
self.assemble_textheader() |
|
199
|
|
|
logging.debug("header:\n%s", self.textheader) |
|
200
|
|
|
|
|
201
|
|
|
# parse geolocation data |
|
202
|
|
|
gr_scia_geo = cl_mds_group.get("geoL_scia") |
|
203
|
|
|
tan_h = gr_scia_geo["tan_h"] |
|
204
|
|
|
# lat and lon are integers in degrees * 10^6 |
|
205
|
|
|
lats_all = gr_scia_geo["tang_ground_point"]["lat"] |
|
206
|
|
|
lons_all = gr_scia_geo["tang_ground_point"]["lon"] |
|
207
|
|
|
sza_all = gr_scia_geo["sun_zen_ang"] |
|
208
|
|
|
saa_all = (gr_scia_geo["sun_azi_ang"] - gr_scia_geo["los_azi_ang"]) |
|
209
|
|
|
sat_h_all = gr_scia_geo["sat_h"] |
|
210
|
|
|
earth_rad_all = gr_scia_geo["earth_rad"] |
|
211
|
|
|
# lat and lon are integers in degrees * 10^6 |
|
212
|
|
|
subsatlat_all = gr_scia_geo["sub_sat_point"]["lat"] |
|
213
|
|
|
subsatlon_all = gr_scia_geo["sub_sat_point"]["lon"] |
|
214
|
|
|
# fix longitudes to [0°, 360°) |
|
215
|
|
|
lons_all[np.where(lons_all < 0)] += 360000000 |
|
216
|
|
|
|
|
217
|
|
|
if cl_n_readouts > 2: |
|
218
|
|
|
tangent_heights = 0.5 * (tan_h[1::cl_n_readouts, 2] + tan_h[2::cl_n_readouts, 0]) |
|
219
|
|
|
tp_lats = 0.5 * (lats_all[1::cl_n_readouts, 2] + lats_all[2::cl_n_readouts, 0]) * 1e-6 |
|
220
|
|
|
tp_lons = 0.5 * (lons_all[1::cl_n_readouts, 2] + lons_all[2::cl_n_readouts, 0]) * 1e-6 |
|
221
|
|
|
sza_toa = 0.5 * (sza_all[1::cl_n_readouts, 2] + sza_all[2::cl_n_readouts, 0]) |
|
222
|
|
|
saa_toa = 0.5 * (saa_all[1::cl_n_readouts, 2] + saa_all[2::cl_n_readouts, 0]) |
|
223
|
|
|
sat_hs = sat_h_all.reshape((-1, cl_n_readouts)).mean(axis=1) |
|
224
|
|
|
earth_rads = earth_rad_all.reshape((-1, cl_n_readouts)).mean(axis=1) |
|
225
|
|
|
subsatlat = subsatlat_all.reshape((-1, cl_n_readouts)).mean(axis=1) * 1e-6 |
|
226
|
|
|
subsatlon = subsatlon_all.reshape((-1, cl_n_readouts)).mean(axis=1) * 1e-6 |
|
227
|
|
|
else: |
|
228
|
|
|
tangent_heights = tan_h[::cl_n_readouts, 1] |
|
229
|
|
|
tp_lats = lats_all[::cl_n_readouts, 1] * 1e-6 |
|
230
|
|
|
tp_lons = lons_all[::cl_n_readouts, 1] * 1e-6 |
|
231
|
|
|
sza_toa = sza_all[::cl_n_readouts, 1] |
|
232
|
|
|
saa_toa = saa_all[::cl_n_readouts, 1] |
|
233
|
|
|
sat_hs = sat_h_all[::cl_n_readouts] |
|
234
|
|
|
earth_rads = earth_rad_all[::cl_n_readouts] |
|
235
|
|
|
subsatlat = subsatlat_all[::cl_n_readouts] * 1e-6 |
|
236
|
|
|
subsatlon = subsatlon_all[::cl_n_readouts] * 1e-6 |
|
237
|
|
|
|
|
238
|
|
|
logging.debug("tangent altitudes: %s", tangent_heights) |
|
239
|
|
|
nalt = len(tangent_heights) |
|
240
|
|
|
|
|
241
|
|
|
centre = _middle_coord(lats_all[0, 1] * 1e-6, lons_all[0, 1] * 1e-6, |
|
242
|
|
|
lats_all[nalt - 2, 1] * 1e-6, lons_all[nalt - 2, 1] * 1e-6) |
|
243
|
|
|
|
|
244
|
|
|
cent_lat_lon = (centre[0], |
|
245
|
|
|
# fix longitudes to [0, 360.) |
|
246
|
|
|
centre[1] if centre[1] >= 0. else 360. + centre[1], |
|
247
|
|
|
lats_all[0, 0] * 1e-6, lons_all[0, 0] * 1e-6, |
|
248
|
|
|
lats_all[0, 2] * 1e-6, lons_all[0, 2] * 1e-6, |
|
249
|
|
|
lats_all[nalt - 2, 0] * 1e-6, lons_all[nalt - 2, 0] * 1e-6, |
|
250
|
|
|
lats_all[nalt - 2, 2] * 1e-6, lons_all[nalt - 2, 2] * 1e-6) |
|
251
|
|
|
|
|
252
|
|
|
toa = 100. |
|
253
|
|
|
# to satellite first |
|
254
|
|
|
los_calc = np.degrees(np.arccos(0.0)) |
|
255
|
|
|
sza_tp_h = sza_toa.copy() |
|
256
|
|
|
saa_tp_h = saa_toa.copy() |
|
257
|
|
|
los_tp_h = np.full_like(tangent_heights, los_calc) |
|
258
|
|
|
los_toa_h = np.full_like(tangent_heights, los_calc) |
|
259
|
|
|
sza_sat_h, los_sat_h, saa_sat_h = _calc_angles( |
|
260
|
|
|
sza_toa, los_calc, saa_toa, |
|
261
|
|
|
tangent_heights, sat_hs, earth_rads) |
|
262
|
|
|
# angles toa |
|
263
|
|
|
los_calc = np.degrees(np.arcsin((tangent_heights + earth_rads) / (toa + earth_rads))) |
|
264
|
|
|
# to tangent point |
|
265
|
|
|
los_toa_l = np.full_like(tangent_heights, los_calc) |
|
266
|
|
|
sza_tp_l, los_tp_l, saa_tp_l = _calc_angles( |
|
267
|
|
|
sza_toa, los_calc, saa_toa, |
|
268
|
|
|
np.full_like(tangent_heights, toa), tangent_heights, |
|
269
|
|
|
earth_rads) |
|
270
|
|
|
# to satellite |
|
271
|
|
|
sza_sat_l, los_sat_l, saa_sat_l = _calc_angles( |
|
272
|
|
|
sza_toa, los_calc, saa_toa, |
|
273
|
|
|
np.full_like(tangent_heights, toa), sat_hs, |
|
274
|
|
|
earth_rads) |
|
275
|
|
|
|
|
276
|
|
|
sza_sat_h[np.where(tangent_heights <= toa)] = 0. |
|
277
|
|
|
sza_sat_l[np.where(tangent_heights > toa)] = 0. |
|
278
|
|
|
saa_sat_h[np.where(tangent_heights <= toa)] = 0. |
|
279
|
|
|
saa_sat_l[np.where(tangent_heights > toa)] = 0. |
|
280
|
|
|
los_sat_h[np.where(tangent_heights <= toa)] = 0. |
|
281
|
|
|
los_sat_l[np.where(tangent_heights > toa)] = 0. |
|
282
|
|
|
|
|
283
|
|
|
sza_tp_h[np.where(tangent_heights <= toa)] = 0. |
|
284
|
|
|
sza_tp_l[np.where(tangent_heights > toa)] = 0. |
|
285
|
|
|
saa_tp_h[np.where(tangent_heights <= toa)] = 0. |
|
286
|
|
|
saa_tp_l[np.where(tangent_heights > toa)] = 0. |
|
287
|
|
|
los_tp_h[np.where(tangent_heights <= toa)] = 0. |
|
288
|
|
|
los_tp_l[np.where(tangent_heights > toa)] = 0. |
|
289
|
|
|
|
|
290
|
|
|
los_toa_h[np.where(tangent_heights <= toa)] = 0. |
|
291
|
|
|
los_toa_l[np.where(tangent_heights > toa)] = 0. |
|
292
|
|
|
|
|
293
|
|
|
sza_sat = sza_sat_h + sza_sat_l |
|
294
|
|
|
saa_sat = saa_sat_h + saa_sat_l |
|
295
|
|
|
los_sat = los_sat_h + los_sat_l |
|
296
|
|
|
|
|
297
|
|
|
sza_tp = sza_tp_h + sza_tp_l |
|
298
|
|
|
saa_tp = saa_tp_h + saa_tp_l |
|
299
|
|
|
los_tp = los_tp_h + los_tp_l |
|
300
|
|
|
|
|
301
|
|
|
los_toa = los_toa_h + los_toa_l |
|
302
|
|
|
|
|
303
|
|
|
logging.debug("TP sza, saa, los: %s, %s, %s", sza_tp, saa_tp, los_tp) |
|
304
|
|
|
logging.debug("TOA sza, saa, los: %s, %s, %s", sza_toa, saa_toa, los_toa) |
|
305
|
|
|
logging.debug("SAT sza, saa, los: %s, %s, %s", sza_sat, saa_sat, los_sat) |
|
306
|
|
|
|
|
307
|
|
|
# save the data to the limb scan class |
|
308
|
|
|
self.nalt = nalt |
|
309
|
|
|
self.orbit_state = (self.metadata["orbit"], state_in_orbit, |
|
310
|
|
|
self.metadata["state_id"], |
|
311
|
|
|
self.metadata["nr_profile"], self.metadata["act_profile"]) |
|
312
|
|
|
self.date = (state_dt.year, state_dt.month, state_dt.day, |
|
313
|
|
|
state_dt.hour, state_dt.minute, state_dt.second) |
|
314
|
|
|
self.orbit_phase = orb_phase |
|
315
|
|
|
self.cent_lat_lon = cent_lat_lon |
|
316
|
|
|
# pre-set the limb_data |
|
317
|
|
|
if self._limb_data_dtype is None: |
|
318
|
|
|
self._limb_data_dtype = _limb_data_dtype[:] |
|
319
|
|
|
self.limb_data = np.zeros((self.nalt), dtype=self._limb_data_dtype) |
|
320
|
|
|
self.limb_data["sub_sat_lat"] = subsatlat |
|
321
|
|
|
self.limb_data["sub_sat_lon"] = subsatlon |
|
322
|
|
|
self.limb_data["tp_lat"] = tp_lats |
|
323
|
|
|
self.limb_data["tp_lon"] = tp_lons |
|
324
|
|
|
self.limb_data["tp_alt"] = tangent_heights |
|
325
|
|
|
self.limb_data["tp_sza"] = sza_tp |
|
326
|
|
|
self.limb_data["tp_saa"] = saa_tp |
|
327
|
|
|
self.limb_data["tp_los"] = los_tp |
|
328
|
|
|
self.limb_data["toa_sza"] = sza_toa |
|
329
|
|
|
self.limb_data["toa_saa"] = saa_toa |
|
330
|
|
|
self.limb_data["toa_los"] = los_toa |
|
331
|
|
|
self.limb_data["sat_sza"] = sza_sat |
|
332
|
|
|
self.limb_data["sat_saa"] = saa_sat |
|
333
|
|
|
self.limb_data["sat_los"] = los_sat |
|
334
|
|
|
self.limb_data["sat_alt"] = sat_hs |
|
335
|
|
|
self.limb_data["earth_rad"] = earth_rads |
|
336
|
|
|
return 0 |
|
337
|
|
|
|
|
338
|
1 |
|
def read_hdf5_limb_state_spectral_data(self, hf, lstate_id, cl_id): |
|
339
|
|
|
"""SCIAMACHY level 1c spectral data |
|
340
|
|
|
|
|
341
|
|
|
Parameters |
|
342
|
|
|
---------- |
|
343
|
|
|
hf : opened file |
|
344
|
|
|
Pointer to the opened level 1c HDF5 file |
|
345
|
|
|
lstate_id : int |
|
346
|
|
|
The limb state id. |
|
347
|
|
|
cl_id : int |
|
348
|
|
|
The spectral cluster number. |
|
349
|
|
|
|
|
350
|
|
|
Returns |
|
351
|
|
|
------- |
|
352
|
|
|
success : int |
|
353
|
|
|
0 on success, |
|
354
|
|
|
1 if an error occurred, for example if the measurement data |
|
355
|
|
|
set for the requested limb and cluster ids is empty. |
|
356
|
|
|
""" |
|
357
|
|
|
cl_mds_group_name = "/MDS/limb_{0:02d}/cluster_{1:02d}".format(lstate_id, cl_id + 1) |
|
358
|
|
|
cl_mds_group = hf.get(cl_mds_group_name) |
|
359
|
|
|
if cl_mds_group is None: |
|
360
|
|
|
return 1 |
|
361
|
|
|
|
|
362
|
|
|
ads_state = hf.get("/ADS/STATES")[lstate_id] |
|
363
|
|
|
cl_n_readouts = ads_state["clus_config"]["num_readouts"][cl_id] |
|
364
|
|
|
|
|
365
|
|
|
pwl = cl_mds_group.get("pixel_wavelength")[:] |
|
366
|
|
|
signal = cl_mds_group.get("pixel_signal")[:] |
|
367
|
|
|
sig_errs = cl_mds_group.get("pixel_signal_error")[:] |
|
368
|
|
|
if cl_n_readouts > 1: |
|
369
|
|
|
# coadd data |
|
370
|
|
|
signal_coadd = signal.reshape((-1, cl_n_readouts, len(pwl))).sum(axis=1) |
|
371
|
|
|
sig_errs = np.sqrt(((sig_errs * signal)**2) |
|
372
|
|
|
.reshape((-1, cl_n_readouts, len(pwl))) |
|
373
|
|
|
.sum(axis=1)) / np.abs(signal_coadd) |
|
374
|
|
|
signal = signal_coadd / cl_n_readouts |
|
375
|
|
|
if np.any(self.wls): |
|
376
|
|
|
# apparently we already have some data, so concatenate |
|
377
|
|
|
self.wls = np.concatenate([self.wls, pwl], axis=0) |
|
378
|
|
|
self._rad_arr = np.concatenate([self._rad_arr, signal], axis=1) |
|
379
|
|
|
self._err_arr = np.concatenate([self._err_arr, sig_errs], axis=1) |
|
380
|
|
|
else: |
|
381
|
|
|
# this seems to be the first time we fill the arrays |
|
382
|
|
|
self.wls = pwl |
|
383
|
|
|
self._rad_arr = signal |
|
384
|
|
|
self._err_arr = sig_errs |
|
385
|
|
|
return 0 |
|
386
|
|
|
|
|
387
|
1 |
|
def read_from_hdf5(self, hf, limb_state_id, state_in_orbit, cluster_ids): |
|
388
|
|
|
"""SCIAMACHY level 1c HDF main interface |
|
389
|
|
|
|
|
390
|
|
|
This should be the function to be used for HDF5 reading. |
|
391
|
|
|
It reads the common data and iterates over the spectral |
|
392
|
|
|
clusters to fill the class data. |
|
393
|
|
|
|
|
394
|
|
|
Parameters |
|
395
|
|
|
---------- |
|
396
|
|
|
hf : opened file |
|
397
|
|
|
Pointer to the opened level 1c HDF5 file |
|
398
|
|
|
limb_state_id : int |
|
399
|
|
|
The limb state id. |
|
400
|
|
|
state_in_orbit : int |
|
401
|
|
|
The number in this batch of states for the header. |
|
402
|
|
|
cluster_ids : int or sequence of ints |
|
403
|
|
|
The spectral cluster numbers. |
|
404
|
|
|
|
|
405
|
|
|
Returns |
|
406
|
|
|
------- |
|
407
|
|
|
success : int |
|
408
|
|
|
0 on success, |
|
409
|
|
|
1 if an error occurred, for example if the measurement data |
|
410
|
|
|
set for the requested limb and cluster ids is empty. |
|
411
|
|
|
""" |
|
412
|
|
|
import numpy.lib.recfunctions as rfn |
|
413
|
|
|
if not hasattr(cluster_ids, '__getitem__'): |
|
414
|
|
|
cluster_ids = [cluster_ids] |
|
415
|
|
|
|
|
416
|
|
|
if self.read_hdf5_limb_state_common_data(hf, limb_state_id, |
|
417
|
|
|
state_in_orbit, cluster_ids[0]): |
|
418
|
|
|
return 1 |
|
419
|
|
|
|
|
420
|
|
|
for cl_id in cluster_ids: |
|
421
|
|
|
self.read_hdf5_limb_state_spectral_data(hf, limb_state_id, cl_id) |
|
422
|
|
|
|
|
423
|
|
|
self.npix = len(self.wls) |
|
424
|
|
|
|
|
425
|
|
|
rads = np.rec.fromarrays([self._rad_arr], |
|
426
|
|
|
dtype=np.dtype([("rad", 'f4', (self.npix,))])) |
|
427
|
|
|
errs = np.rec.fromarrays([self._err_arr], |
|
428
|
|
|
dtype=np.dtype([("err", 'f4', (self.npix,))])) |
|
429
|
|
|
self.limb_data = rfn.merge_arrays([self.limb_data, rads, errs], |
|
430
|
|
|
usemask=False, asrecarray=True, flatten=True) |
|
431
|
|
|
self._limb_data_dtype = self.limb_data.dtype |
|
432
|
|
|
|
|
433
|
|
|
return 0 |
|
434
|
|
|
|
st-bender /
sciapy
