-
Notifications
You must be signed in to change notification settings - Fork 1
/
sediment_indices.py
155 lines (121 loc) · 4.93 KB
/
sediment_indices.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
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
# Pull images and determine sediment indices for AOIs
# Remember to update ee.Initialize() with authenticated project - LINE 18
# Initial attempt in GEE:
# https://code.earthengine.google.com/77a2e41436a615fe1c04c21fec61e9d0
# https://developers.google.com/earth-engine/guides/python_install
# pip install earthengine-api --upgrade
# https://developers.google.com/earth-engine/guides/auth
import ee
import numpy as np
import matplotlib.pyplot as plt
import shapefile
import pandas as pd
#ee.Authenticate(force=True)
ee.Authenticate()
ee.Initialize(project = 'ee-fortschthomas52')
#ee.Initialize(project = 'ee-davidmkahler-limpopo')
# Define analysis area
# the analysis area is transformed into a rectangle in the .sampleRectangle/.getInfo phases
# enter shapefiles here:
#sf = shapefile.Reader("/Volumes/dmk/gis/limpopo/kruger/logger_sites/reference_polygons/balule/balule.shp")
# sf = shapefile.Reader("/Users/hydro3/Documents/KrugerSensors/reference_polygons/balule/balule.shp")
# shapes = sf.shapes()
# points = shapes[0].points
# aoi = ee.Geometry.Polygon(points)
aoi = ee.Geometry.Polygon(
[[31.717329298139624, -24.055719459004635],
[31.716846500516944, -24.057384943478993],
[31.718235884786658, -24.058736909100467],
[31.72093418705564, -24.05800214694256],
[31.71991494762998, -24.054984679571938]], None, False)
start = pd.date_range(start= '2023-04-01' , end='2024-05-17' ,
freq='5d')
end = pd.date_range(start='2023-04-06' , end='2024-05-22' ,
freq='5d')
dates = pd.DataFrame ({'start': start , 'end': end})
def pulldata(startDate, endDate):
# Define source data
image = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')\
.filterDate(startDate, endDate)\
.select('B2', 'B3', 'B4', 'B8', 'B11')
#.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 10)\
# CRS is not the same.
# proj = image.first().select('B2').projection() # EPSG:32656, UTM zone 56N (Siberia?)
# proj = balule.projection() # EPSG:4326
# Export arrays
# https://gist.github.com/jdbcode/f4d56d72f7fc5beeaa3859999b1f5c3d
# https://gist.github.com/jdbcode/f4d56d72f7fc5beeaa3859999b1f5c3d?permalink_comment_id=3355627#gistcomment-3355627
mosaic = image.median().reproject(crs='EPSG:32736', scale=1) # This allows us to set the resolution.
band_arrs = mosaic.sampleRectangle(region=aoi)
# Get individual band arrays.
band_arr_b2 = band_arrs.get('B2') # Blue
band_arr_b3 = band_arrs.get('B3') # Green
band_arr_b4 = band_arrs.get('B4') # Red
band_arr_b8 = band_arrs.get('B8') # NIR
band_arr_b11 = band_arrs.get('B11') # SWIR
# Transfer the arrays from server to client and cast as np array.
b2 = np.array(band_arr_b2.getInfo()) # b2 Blue
b3 = np.array(band_arr_b3.getInfo()) # b3 Green
b4 = np.array(band_arr_b4.getInfo()) # b4 Red
b8 = np.array(band_arr_b8.getInfo()) # b8 NIR
b11 = np.array(band_arr_b11.getInfo()) # b11 SWIR
np_arr_b4 = np.expand_dims(b4, 2)
np_arr_b3 = np.expand_dims(b3, 2)
np_arr_b2 = np.expand_dims(b2, 2)
rgb_img = np.concatenate((np_arr_b4, np_arr_b3, np_arr_b2), 2)
rgb_img = (255*((rgb_img)/3000)).astype('uint8')
plt.imshow(rgb_img)
plt.show()
# Normalized Difference Water Index (NDWI)
# NDWI = ( G - NIR ) / ( G + NIR )
ndwiG = (b3-b8)/(b3+b8) # Gao
# NDWI = ( NIR - SWIR ) / ( NIR + SWIR )
ndwiM = (b8-b11)/(b8+b11) # McFeeters
# NDWI = ( G - SWIR ) / ( G + SWIR )
mndwi = (b3-b11)/(b3+b11) # Modified NDWI
water = ndwiG > -0.02
TSS1 = np.NAN
Secchi1 = np.NAN
TSS2 = np.NAN
Secchi2 = np.NAN
TSS3 = np.NAN
Secchi3 = np.NAN
TSS4 = np.NAN
Ratio = np.NAN
TSS1 = (b3+b4)/2
TSS1 = TSS1 * water
#plt.imshow(TSS1 , cmap= "summer")
#plt.colorbar(orientation="vertical")
#plt.show
TSS1 = np.sum(TSS1) / np.sum(TSS1>0)
#print(TSS1)
Secchi1 = (b2/b4)
Secchi1 = Secchi1 * water
Secchi1 = np.sum(Secchi1)/ np.sum(Secchi1>0)
TSS2 = (b3/b4)
TSS2 = TSS2 * water
TSS2 = np.sum(TSS2) / np.sum(TSS2>0)
Secchi2 = (b4/b3)
Secchi2 = Secchi2 * water
Secchi2 = np.sum(Secchi2)/ np.sum(Secchi2>0)
TSS3 = (b8/b3 , b8/b4)
TSS3 = TSS3 * water
TSS3 = np.sum(TSS3) / np.sum(TSS3>0)
Secchi3 = (b4/b2)+b2
Secchi3 = Secchi3 * water
Secchi3 = np.sum(Secchi3)/ np.sum(Secchi3>0)
TSS4 = (b4/b3)+b8
TSS4 = TSS4 * water
TSS4 = np.sum(TSS4) / np.sum(TSS4>0)
Ratio = (ndwiG/ndwiM)
Ratio = Ratio * water
Ratio = np.sum(Ratio)/ np.sum(Ratio>0)
f = open("Sediment_Indices5dayTest.txt", "a")
f.write(str(startDate) + ", " + str(endDate) + ", " + str(TSS1) + ", " +
str(Secchi1) + ", " +
str(TSS2) + ", " + str(Secchi2) + ", " +
str(TSS3) + ", " + str(Secchi3) + ", " +
str(TSS4) + "," + str(Ratio) + '\n')
f.close()
for i in range(len(dates)):
pulldata(dates["start"][i], dates["end"][i])