-
Notifications
You must be signed in to change notification settings - Fork 6
/
rms_indexes.py
102 lines (87 loc) · 4.01 KB
/
rms_indexes.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
import ee
from zipfile import ZipFile
from io import BytesIO
import os
import requests
class S2indexes:
def __init__(self, area, dir, date_from, date_end, scope):
"""
given an area defined by a geoJSON, it returns rasters of
remote sensing indexes at the specified date at granularuity defined by the scope parameter
Args:
area: geoJSON, use squaretogeojson to generate
dir: directory where to save the easter
date_from, date_end: nightlights for what point in time?
scope (str): country or urban?
"""
self.area = area
self.dir = dir
self.date_from = date_from
self.date_end = date_end
self.scope = scope
self.files = None
def download(self):
print('INFO: downloading rms indexes for area of interest ...')
if os.path.exists(self.dir + str(self.area["coordinates"]) + "NDVI_max.tif") \
and os.path.exists(self.dir + str(self.area["coordinates"]) + "NDBI_max.tif") \
and os.path.exists(self.dir + str(self.area["coordinates"]) + "NDWI_max.tif"):
self.files = [str(self.area["coordinates"]) + "NDVI_max.tif", str(self.area["coordinates"]) + "NDBI_max.tif", str(self.area["coordinates"]) + "NDWI_max.tif"]
print('INFO: NDs data for {} already downloaded'.format(self.area["coordinates"]))
else:
ee.Initialize()
GREEN = 'B3'
RED = 'B4'
NIR = 'B8'
SWIR = 'B11'
sentinel = ee.ImageCollection('COPERNICUS/S2') \
.filterDate(ee.Date(self.date_from), ee.Date(self.date_end)) \
.filterBounds(self.area) \
.select([GREEN, RED, NIR, SWIR]) \
.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 70)
def addIndices(image):
ndvi = image.normalizedDifference([NIR, RED])
ndbi = image.normalizedDifference([SWIR, NIR])
ndwi = image.normalizedDifference([GREEN, NIR])
return image.addBands(ndvi.rename('NDVI')) \
.addBands(ndbi.rename('NDBI')) \
.addBands(ndwi.rename('NDWI'))
img = sentinel.map(addIndices).select(['NDVI', 'NDBI', 'NDWI']).reduce("max").clip(self.area)
# download
if self.scope == 'urban':
print('INFO: NDs scope > urban')
scale = 100
else:
print('INFO: NDs scope -> country')
scale = 5000
for b in ['NDVI_max', 'NDBI_max', 'NDWI_max']:
url = img.select(b).getDownloadUrl({'crs': 'EPSG:4326', 'region': self.area, 'scale': scale})
r = requests.get(url)
z = ZipFile(BytesIO(r.content))
z.extract(z.namelist()[1], self.dir)
os.rename(self.dir + z.namelist()[1], self.dir + str(self.area["coordinates"]) + b+'.tif')
self.files = [str(self.area["coordinates"]) + "NDVI_max.tif", str(self.area["coordinates"]) + "NDBI_max.tif",
str(self.area["coordinates"]) + "NDWI_max.tif"]
def rms_values(self, longitudes, latitudes):
"""
Given a dataset with latitude and longitude columns, it returns the nightlight value at each point.
Args:
longitudes: list of longitudes
latitudes: list of latitudes
Returns:
Series
"""
import rasterio
try:
NDVI = rasterio.open(self.dir + self.files[0])
NDBI = rasterio.open(self.dir + self.files[1])
NDWI = rasterio.open(self.dir + self.files[2])
except MemoryError:
print('Remote sensing indexes rasters too big!')
raise
veg, build, wat = [], [], []
for lon, lat in zip(longitudes, latitudes):
i, j = NDVI.index(lon, lat)
veg.append(NDVI.read(1)[i, j])
build.append(NDBI.read(1)[i, j])
wat.append(NDWI.read(1)[i, j])
return veg, build, wat