|
1
|
|
|
#!/usr/bin/env python |
|
2
|
|
|
# -*- coding: utf-8 -*- |
|
3
|
|
|
|
|
4
|
|
|
# Copyright (c) 2017 Adam.Dybbroe |
|
5
|
|
|
|
|
6
|
|
|
# Author(s): |
|
7
|
|
|
|
|
8
|
|
|
# Adam.Dybbroe <[email protected]> |
|
9
|
|
|
# Thomas Lepellt, DWD |
|
10
|
|
|
|
|
11
|
|
|
# This program is free software: you can redistribute it and/or modify |
|
12
|
|
|
# it under the terms of the GNU General Public License as published by |
|
13
|
|
|
# the Free Software Foundation, either version 3 of the License, or |
|
14
|
|
|
# (at your option) any later version. |
|
15
|
|
|
|
|
16
|
|
|
# This program is distributed in the hope that it will be useful, |
|
17
|
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
18
|
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
19
|
|
|
# GNU General Public License for more details. |
|
20
|
|
|
|
|
21
|
|
|
# You should have received a copy of the GNU General Public License |
|
22
|
|
|
# along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
23
|
|
|
|
|
24
|
|
|
"""Depending on the satellite sensor spectral band the upwelling infrared |
|
25
|
|
|
radiation in a cloudfree atmosphere will be subject to absorption by |
|
26
|
|
|
atmospherical molecules, mainly water vapour, Ozone, CO2 and other trace |
|
27
|
|
|
gases. This absorption is depending on the observational path length. Usually |
|
28
|
|
|
the satellite signal is subject to what is referred to as limb cooling, such |
|
29
|
|
|
that observed brightness temperatures under large view angles are lower than |
|
30
|
|
|
those under smaller view angles. |
|
31
|
|
|
|
|
32
|
|
|
This code seek to correct for this cooling, so that it is possible to get a |
|
33
|
|
|
more consistent looking satellite image also when using imager bands which are |
|
34
|
|
|
not true window channels, but rather subject to some absorption ('dirty window |
|
35
|
|
|
channel'). |
|
36
|
|
|
|
|
37
|
|
|
The goal is eventually to do this atmospheric correction with the use of off |
|
38
|
|
|
line Radiative Transfer Model simulations giving look up tables of atmospheric |
|
39
|
|
|
absorption for various standard atmospheres and satellite zenith angles and |
|
40
|
|
|
spectral bands. |
|
41
|
|
|
|
|
42
|
|
|
But, for now, we apply and old parametric method implemented at DWD since long |
|
43
|
|
|
ago: |
|
44
|
|
|
|
|
45
|
|
|
J.Asmus (August 2017): 'Some information to our atmospheric correction. This |
|
46
|
|
|
routine is very old, as Katja told. It running in DWD since around 1989 at |
|
47
|
|
|
first on DEC Micro VAX 4000-300, later on SGI workstation (O2) both as FORTRAN |
|
48
|
|
|
program, later in CineSat and now in PyTroll. This correction is very simple |
|
49
|
|
|
but it had to run on a Micro VAX. The idea was to get some more realistic |
|
50
|
|
|
temperatures in the higher latitudes. The source for this atmospheric |
|
51
|
|
|
correction must be a paper from the US Air Force in the 1960 or 1970 years |
|
52
|
|
|
(?). Unfortunaly I have no more information about the source and the colleagues |
|
53
|
|
|
from this time are not longer at DWD.' |
|
54
|
|
|
|
|
55
|
|
|
""" |
|
56
|
|
|
|
|
57
|
|
|
|
|
58
|
|
|
import numpy as np |
|
59
|
|
|
import logging |
|
60
|
|
|
|
|
61
|
|
|
LOG = logging.getLogger(__name__) |
|
62
|
|
|
|
|
63
|
|
|
|
|
64
|
|
|
class AtmosphericalCorrection(object): |
|
65
|
|
|
|
|
66
|
|
|
"""Container for the IR atmospherical correction of satellite imager IR (dirty) |
|
67
|
|
|
window bands |
|
68
|
|
|
|
|
69
|
|
|
""" |
|
70
|
|
|
|
|
71
|
|
|
def __init__(self, platform_name, sensor, **kwargs): |
|
72
|
|
|
self.platform_name = platform_name |
|
73
|
|
|
self.sensor = sensor |
|
74
|
|
|
self.coeff_filename = None |
|
75
|
|
|
|
|
76
|
|
|
if 'atmosphere' in kwargs: |
|
77
|
|
|
atm_type = kwargs['atmosphere'] |
|
78
|
|
|
else: |
|
79
|
|
|
atm_type = None |
|
80
|
|
|
|
|
81
|
|
|
if atm_type: |
|
82
|
|
|
raise AttributeError("Atmosphere type %s not supported", atm_type) |
|
83
|
|
|
|
|
84
|
|
|
LOG.info("Atmospherical correction by old DWD parametric method...") |
|
85
|
|
|
|
|
86
|
|
|
def get_correction(self, sat_zenith, bandname, data): |
|
87
|
|
|
""" |
|
88
|
|
|
Get the correction depending on satellite zenith angle and band |
|
89
|
|
|
""" |
|
90
|
|
|
|
|
91
|
|
|
# From band name we will eventually need the effective wavelength |
|
92
|
|
|
# and then use that to find the atm contribution depedning on zenith |
|
93
|
|
|
# angle from a LUT |
|
94
|
|
|
|
|
95
|
|
|
return viewzen_corr(data, sat_zenith) |
|
96
|
|
|
|
|
97
|
|
|
|
|
98
|
|
|
def viewzen_corr(data, view_zen): |
|
99
|
|
|
"""Apply atmospheric correction on the given *data* using the |
|
100
|
|
|
specified satellite zenith angles (*view_zen*). Both input data |
|
101
|
|
|
are given as 2-dimensional Numpy (masked) arrays, and they should |
|
102
|
|
|
have equal shapes. |
|
103
|
|
|
The *data* array will be changed in place and has to be copied before. |
|
104
|
|
|
""" |
|
105
|
|
|
def ratio(value, v_null, v_ref): |
|
106
|
|
|
return (value - v_null) / (v_ref - v_null) |
|
107
|
|
|
|
|
108
|
|
|
def tau0(t): |
|
109
|
|
|
T_0 = 210.0 |
|
110
|
|
|
T_REF = 320.0 |
|
111
|
|
|
TAU_REF = 9.85 |
|
112
|
|
|
return (1 + TAU_REF)**ratio(t, T_0, T_REF) - 1 |
|
113
|
|
|
|
|
114
|
|
|
def tau(t): |
|
115
|
|
|
T_0 = 170.0 |
|
116
|
|
|
T_REF = 295.0 |
|
117
|
|
|
TAU_REF = 1.0 |
|
118
|
|
|
M = 4 |
|
119
|
|
|
return TAU_REF * ratio(t, T_0, T_REF)**M |
|
120
|
|
|
|
|
121
|
|
|
def delta(z): |
|
122
|
|
|
Z_0 = 0.0 |
|
123
|
|
|
Z_REF = 70.0 |
|
124
|
|
|
DELTA_REF = 6.2 |
|
125
|
|
|
return (1 + DELTA_REF)**ratio(z, Z_0, Z_REF) - 1 |
|
126
|
|
|
|
|
127
|
|
|
y0, x0 = np.ma.where(view_zen == 0) |
|
128
|
|
|
data[y0, x0] += tau0(data[y0, x0]) |
|
129
|
|
|
|
|
130
|
|
|
y, x = np.ma.where((view_zen > 0) & (view_zen < 90) & (~data.mask)) |
|
131
|
|
|
data[y, x] += tau(data[y, x]) * delta(view_zen[y, x]) |
|
132
|
|
|
return data |
|
133
|
|
|
|
|
134
|
|
|
|
|
135
|
|
|
if __name__ == "__main__": |
|
136
|
|
|
|
|
137
|
|
|
this = AtmosphericalCorrection('Suomi-NPP', 'viirs') |
|
138
|
|
|
SHAPE = (1000, 3000) |
|
139
|
|
|
NDIM = SHAPE[0] * SHAPE[1] |
|
140
|
|
|
SATZ = np.ma.arange(NDIM).reshape(SHAPE) * 60. / float(NDIM) |
|
141
|
|
|
TBS = np.ma.arange(NDIM).reshape(SHAPE) * 80.0 / float(NDIM) + 220. |
|
142
|
|
|
atm_corr = this.get_correction(SATZ, 'M4', TBS) |
|
143
|
|
|
|
pytroll /
pyspectral
