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
|
#!/usr/bin/env python3
# This program is a direct translation of the sample program
# "write_ndvi.c" bundled with the EPR-API distribution.
#
# Source code of the C program is available at:
# https://github.com/bcdev/epr-api/blob/master/src/examples/write_ndvi.c
"""Example for using the epr-api.
Demonstrates how to open a MERIS L1b product and calculate the NDVI.
This example does not demonstrate how to write good and safe code.
It is reduced to the essentials for working with the epr-api.
Calling sequence::
$ python3 write_ndvi.py <envisat-product> <output-file>
for example::
$ python3 write_ndvi.py MER_RR__1P_test.N1 my_ndvi.raw
"""
import sys
import struct
import logging
import epr
def main(*argv):
if not argv:
argv = sys.argv
if len(argv) != 3:
print(
""""
Usage:
python3 write_ndvi <envisat-product> <output-file>
where envisat-product is the input filename
and output-file is the output filename.
Example:
MER_RR__1P_TEST.N1 my_ndvi.raw
"""
)
sys.exit(1)
# Open the product
with epr.open(argv[1]) as product:
# The NDVI shall be calculated using bands 6 and 8.
band1_name = "radiance_6"
band2_name = "radiance_10"
band1 = product.get_band(band1_name)
band2 = product.get_band(band2_name)
# Allocate memory for the rasters
width = product.get_scene_width()
height = product.get_scene_height()
subsampling_x = 1
subsampling_y = 1
raster1 = band1.create_compatible_raster(
width, height, subsampling_x, subsampling_y
)
raster2 = band2.create_compatible_raster(
width, height, subsampling_x, subsampling_y
)
# Read the radiance into the raster.
offset_x = 0
offset_y = 0
logging.info("read %r data", band1_name)
band1.read_raster(offset_x, offset_y, raster1)
logging.info("read %r data", band2_name)
band2.read_raster(offset_x, offset_y, raster2)
# Open the output file
logging.info("write ndvi to %r", argv[2])
with open(argv[2], "wb") as out_stream:
# Loop over all pixel and calculate the NDVI.
#
# @NOTE: looping over data matrices is not the best solution.
# It is done here just for demostrative purposes
for j in range(height):
for i in range(width):
rad1 = raster1.get_pixel(i, j)
rad2 = raster2.get_pixel(i, j)
if (rad1 + rad2) != 0.0:
ndvi = (rad2 - rad1) / (rad2 + rad1)
else:
ndvi = -1.0
out_stream.write(struct.pack("f", ndvi))
logging.info("ndvi was written success")
if __name__ == "__main__":
main()
|
