This repository has been archived by the owner on Jun 14, 2021. It is now read-only.
/
l8_ndvi.py
109 lines (78 loc) · 3.8 KB
/
l8_ndvi.py
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
import json
import base64
from io import BytesIO
import numpy as np
from PIL import Image
import rasterio as rio
from rasterio import warp
from rasterio.enums import Resampling
from rio_toa import reflectance
from remotepixel import utils
def point(scene, coord):
scene_params = utils.landsat_parse_scene_id(scene)
meta_data = utils.landsat_get_mtl(scene)
landsat_address = f's3://landsat-pds/{scene_params["key"]}'
E = float(utils.landsat_mtl_extract(meta_data, 'SUN_ELEVATION'))
MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_4'))
AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_4'))
band_address = f'{landsat_address}_B4.TIF'
with rio.open(band_address) as band:
lon_srs, lat_srs = warp.transform('EPSG:4326', band.crs, [coord[0]], [coord[1]])
b4 = list(band.sample([(lon_srs[0], lat_srs[0])]))[0]
b4 = reflectance.reflectance(b4, MR, AR, E, src_nodata=0)[0]
MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_5'))
AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_5'))
band_address = f'{landsat_address}_B5.TIF'
with rio.open(band_address) as band:
lon_srs, lat_srs = warp.transform('EPSG:4326', band.crs, [coord[0]], [coord[1]])
b5 = list(band.sample([(lon_srs[0], lat_srs[0])]))[0]
b5 = reflectance.reflectance(b5, MR, AR, E, src_nodata=0)[0]
if b4 * b5 > 0:
ratio = (b5 - b4) / (b5 + b4)
else:
ratio = 0
out = {
'ndvi': ratio,
'date': scene_params['date'],
'cloud': utils.landsat_mtl_extract(meta_data, 'CLOUD_COVER')
}
return out
def area(scene, bbox):
max_width = 512
max_height = 512
scene_params = utils.landsat_parse_scene_id(scene)
meta_data = utils.landsat_get_mtl(scene)
landsat_address = f's3://landsat-pds/{scene_params["key"]}'
E = float(utils.landsat_mtl_extract(meta_data, 'SUN_ELEVATION'))
MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_4'))
AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_4'))
band_address = f'{landsat_address}_B4.TIF'
with rio.open(band_address) as band:
crs_bounds = warp.transform_bounds('EPSG:4326', band.crs, *bbox)
window = band.window(*crs_bounds, boundless=True)
width = window.num_cols if window.num_cols < max_width else max_width
height = window.num_rows if window.num_rows < max_width else max_width
b4 = band.read(window=window,
out_shape=(height, width), indexes=1,
resampling=Resampling.bilinear, boundless=True)
b4 = reflectance.reflectance(b4, MR, AR, E, src_nodata=0)
MR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_MULT_BAND_5'))
AR = float(utils.landsat_mtl_extract(meta_data, 'REFLECTANCE_ADD_BAND_5'))
band_address = f'{landsat_address}_B5.TIF'
with rio.open(band_address) as band:
crs_bounds = warp.transform_bounds('EPSG:4326', band.crs, *bbox)
window = band.window(*crs_bounds, boundless=True)
width = window.num_cols if window.num_cols < max_width else max_width
height = window.num_rows if window.num_rows < max_width else max_width
b5 = band.read(window=window,
out_shape=(height, width), indexes=1,
resampling=Resampling.bilinear, boundless=True)
b5 = reflectance.reflectance(b5, MR, AR, E, src_nodata=0)
ratio = np.where( (b5 * b4) > 0, np.nan_to_num((b5 - b4) / (b5 + b4)), -1)
ratio = np.where(ratio > -1,
utils.linear_rescale(ratio, in_range=[-1,1], out_range=[0, 255]), 0)
img = Image.fromarray(ratio).convert('RGB')
sio = BytesIO()
img.save(sio, 'jpeg', subsampling=0, quality=100)
sio.seek(0)
return base64.b64encode(sio.getvalue()).decode()