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
|
NDVI computation
----------------
This tutorial shows how to use PyEPR_ to open a MERIS_ L1B product, compute
the *Normalized Difference Vegetation Index* (NDVI) and store it into a flat
binary file.
The example code (:download:`examples/write_ndvi.py`) is a direct
translation of the C sample program `write_ndvi.c`_ bundled with the
EPR API distribution.
The program is invoked as follows:
.. code-block:: sh
$ python write_ndvi.py <envisat-oroduct> <output-file>
.. _PyEPR: https://github.com/avalentino/pyepr
.. _MERIS: http://envisat.esa.int/handbooks/meris
.. _`write_ndvi.c`: https://github.com/bcdev/epr-api/blob/master/src/examples/write_ndvi.c
The code have been kept very simple and it consists in a single function
(:func:`main`) that also performs a minimal command line arguments handling.
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:linenos:
The ENVISAT_ :class:`epr.Product` is opened using the :func:`epr.open`
function.
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 49
The name of the product is in the first argument of the program call.
In order to keep the code simple no check is performed to ensure that the
product is a valid L1B product.
The NDVI is calculated using bands 6 and 8 (the names of these bands are
"radiance_6" and "radiance_10").
:class:`epr.Band` objects are retrieved using the :meth:`epr.Product.get_band`
method:
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 51-56
*band1* and *band2* are used to read the calibrated radiances into the
:class:`epr.Raster` objects that allow to access data matrices with the
radiance values.
Before reading data into the :class:`epr.Raster` objects they have to be
instantiated specifying their size and data type in order to allow the library
to allocate the correct amount of memory.
For sake of simplicity :class:`epr.Raster` object are created with the same
size of the whole product (with no sub-sampling) using the
:meth:`epr.Band.create_compatible_raster` method of the :class:`epr.Band`
class:
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 58-66
Then data are actually loaded into memory using the
:meth:`epr.Band.read_raster` method.
Since :class:`epr.Raster` objects have been defined to match the whole
product, offset parameters are set to zero (data are read starting from
specified offset):
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 68-76
.. note::
in this simplified example it is assumed that there is enough system
memory to hold the two :class:`epr.Raster` objects.
After opening (in binary mode) the stream for the output
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 78-80
the program simply loops over all pixel and calculate the NDVI with the
following formula:
.. math:: NDVI = \frac{radiance_{10} - radiance_8}{radiance_{10} + radiance_8}
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 82-95
This part of the code tries to mimic closely the original C code
(`write_ndvi.c`_)
.. code-block:: c
out_stream = fopen(argv[2], "wb");
for (j = 0; j < height; ++j) {
for (i = 0; i < width; ++i) {
rad1 = epr_get_pixel_as_float(raster1, i, j);
rad2 = epr_get_pixel_as_float(raster2, i, j);
if ((rad1 + rad2) != 0.0) {
ndvi = (rad2 - rad1) / (rad2 + rad1);
} else {
ndvi = -1.0;
}
status = fwrite( & ndvi, sizeof(float), 1, out_stream);
}
}
epr_log_message(e_log_info, "ndvi was written success");
and uses the :meth:`epr.Raster.get_pixel` method to access pixel values and
perform computation.
The Python_ :func:`struct.pack` function together with :meth:`file.write` is
used to write the NDVI of the pixel n the file in binary format.
.. raw:: latex
\fvset{fontsize=\footnotesize}
.. literalinclude:: examples/write_ndvi.py
:lines: 94
.. note::
the entire solution is quite not pythonic_. As an alternative
implementation it could be used the :class:`numpy.ndarray` interface of
:class:`epr.Raster` objects available via the :data:`epr.Raster.data`
property. The NDVI index is computed on all pixels altogether using
vectorized expressions::
# Initialize the entire matrix to -1
ndvi = numpy.zeros((height, width), 'float32') - 1
aux = raster2.data + raster1.data
# indexes of pixel with non null denominator
idx = numpy.where(aux != 0)
# actual NDVI computation
ndvi[idx] = (raster2.data[idx] - raster1.data[idx]) / aux[idx]
Finally data can be saved to file simply using the
:meth:`numpy.ndarray.tofile` method::
ndvi.tofile(out_stream)
.. _ENVISAT: http://envisat.esa.int
.. _Python: http://www.python.org
.. _pythonic: http://www.cafepy.com/article/be_pythonic
|
