-
Notifications
You must be signed in to change notification settings - Fork 1
/
odanah
134 lines (105 loc) · 4.2 KB
/
odanah
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
import ee
import numpy as np
import matplotlib.pyplot as plt
import shapefile
import pandas as pd
#ee.Authenticate(force=True)
# Token generated with all permissions 29 May 2024, DMK
ee.Authenticate()
ee.Initialize(project = 'ee-fortschthomas52')
sf = shapefile.Reader("/Users/hydro3/Documents/USGSSites/shape_files/wisconsin/odanah.shp")
shapes = sf.shapes()
points = shapes[0].points
aoi = ee.Geometry.Polygon(points)
start = pd.date_range(start= '2017-10-01' , end='2024-06-22' ,
freq='5d')
end = pd.date_range(start='2017-10-05' , end='2024-06-26' ,
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=10) # This allows us to set the resolution.
band_arrs = mosaic.sampleRectangle(region=aoi, defaultValue = 0)
# 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
if(np.max(b2)!=0):
#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("Odanah5day.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])