An error occurred while loading the file. Please try again.
-
Thibault Hallouin authoredUltimately, the objective is for the user to be able to get the raw sampled metric values, or the mean and standard deviation of the sampled metric values, or a series of quantiles of the sampled metric values. There are still problems with the standard deviation on rtensor, and the computation of the quantiles does not work on n-dim expressions yet. So the second and third options are not possible yet, so only the raw values can be returned. Nonetheless, the machinery and the choice of where to introduce the summary functionality could be implemented, which is the purpose of this commit. A new parameter of the bootstrap experiment called "summary" is added: it can be given a value of 0 (to get the raw values). In the future, it would also take a value of 1 for mean+std, and 2 for quantiles.
d37e69d2
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
import datetime
import getopt
import os
import platform
import stat
import subprocess
import sys
import glob
import ogr
import osr
from getSentinel2 import readS2TileInfo
from mtdUtils import setNoDataValue,randomword
def findtiles(shp):
s2tg = ogr.Open('./aux_data/s2_tiling_grid.shp')
clip = ogr.Open(shp)
s2tg_ly = s2tg.GetLayer(0)
clip_ly = clip.GetLayer(0)
tiles = []
s2tg_srs = s2tg_ly.GetSpatialRef()
clip_srs = clip_ly.GetSpatialRef()
tr = osr.CoordinateTransformation(clip_srs, s2tg_srs)
cover = None
covl = [0]
for fr in s2tg_ly:
clip_ly.ResetReading()
gr = fr.GetGeometryRef()
for f in clip_ly:
g = f.GetGeometryRef()
g.Transform(tr)
if g.Intersects(gr):
intsect = g.Intersection(gr)
if cover is None:
cover = intsect
else:
cover = cover.Union(intsect)
covl.append(cover.GetArea() / g.GetArea())
tl = fr.GetFieldAsString('Name')
if tl not in tiles:
tiles.append(tl)
dovl = [covl[i] - covl[i - 1] for i in range(1, len(covl))]
idx = sorted(range(len(dovl)), key=lambda k: dovl[k], reverse=True)
tot = 0.0
i = 0
while i < len(dovl) and tot < 1.0:
tot += dovl[idx[i]]
i += 1
final_tiles = [tiles[k] for k in idx[0:i]]
return final_tiles
def main(argv):
try:
opts, args = getopt.getopt(argv, 'o:crm:', ['min-coverage=', 'max-clouds=', 'scan-footprint'])
except getopt.GetoptError as err:
sys.exit(err)
mindc = 0.0
maxcc = 100.0
clipping = False
scan_foot = False
# tile
if os.path.exists(args[0]):
tiles = findtiles(args[0])
else:
tiles = args[0].split(':')
# Parse options
opt_str = ''
for opt, val in opts:
if opt == '-o':
opt_str = opt_str + ' -o \"' + val + '\"'
elif opt == '-c' and os.path.exists(args[0]):
clipping = True
opt_str = opt_str + ' -c \"' + args[0] + '\"'
elif opt == '-r':
opt_str = opt_str + ' -r'
elif opt == '-m':
opt_str = opt_str + ' -m ' + val
elif opt == '--min-coverage':
mindc = float(val)
elif opt == '--max-clouds':
maxcc = float(val)
elif opt == '--scan-footprint':
scan_foot = True
if platform.system() == 'Windows':
f = open('S2DownloadScript.bat', 'w')
f.write("@echo off\n")
polycmd = 'gdal_polygonize '
else:
f = open('S2DownloadScript.sh', 'w')
polycmd = 'gdal_polygonize.py '
g = open('S2SearchLog.txt', 'w')
DEVN = open(os.devnull, 'w')
for tile in tiles:
mrd = tile[0:2]
tl1 = tile[2]
tl2 = tile[3:5]
# date and date-range
sd = datetime.date(int(args[1][0:4]), int(args[1][4:6]), int(args[1][6:8]))
ed = datetime.date(int(args[2][0:4]), int(args[2][4:6]), int(args[2][6:8])) + datetime.timedelta(1)
cd = sd
valid_dates = []
print 'Scanning date interval for tile ' + tile + '...'
while cd != ed:
yy = str(cd.year)
mm = str(cd.month)
dd = str(cd.day)
page = '/'.join([mrd, tl1, tl2, yy, mm, dd, '0'])
url = 'http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/' + page + '/tileInfo.json'
ret = subprocess.call('wget -q ' + url, shell=True)
if ret == 0:
# Patch: read over all possible "sequences" of a same date
good_seq = 0
good_info = readS2TileInfo('tileInfo.json')
os.remove('tileInfo.json')
curr_seq = 1
while ret == 0:
page = page[:-1] + str(curr_seq)
url = 'http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/' + page + '/tileInfo.json'
ret = subprocess.call('wget -q ' + url, shell=True)
if ret == 0:
info = readS2TileInfo('tileInfo.json')
os.remove('tileInfo.json')
if info['dataCoveragePercentage'] > good_info['dataCoveragePercentage']:
good_info = info
good_seq = curr_seq
curr_seq += 1
info = good_info
datestr = ('%04d%02d%02d' % (cd.year, cd.month, cd.day), good_seq)
if info['dataCoveragePercentage'] >= mindc and info['cloudyPixelPercentage'] <= maxcc:
# Check image footprint versus shape
clip_ok = True
if clipping and scan_foot:
fn = randomword(16)
page = page[:-1] + str(good_seq)
url = 'http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/tiles/' + page + '/B01.jp2'
subprocess.call('wget -q ' + url, shell=True)
subprocess.call('otbcli_BandMath -il B01.jp2 -exp \"im1b1>0\" -out ' + fn + '.tif uint8',
shell=True, stdout=DEVN)
setNoDataValue(fn + '.tif')
subprocess.call(polycmd + fn + '.tif -f \"ESRI Shapefile" ' + fn + '.shp',
shell=True, stdout=DEVN)
ds_query = ogr.Open(args[0])
ds_footp = ogr.Open(fn + '.shp')
ly_query = ds_query.GetLayer(0)
ly_footp = ds_footp.GetLayer(0)
ly_query_srs = ly_query.GetSpatialRef()
ly_footp_srs = ly_footp.GetSpatialRef()
ct = osr.CoordinateTransformation(ly_footp_srs, ly_query_srs)
intsec = False
for fq in ly_query:
for ff in ly_footp:
gf = ff.GetGeometryRef()
gf.Transform(ct)
if fq.GetGeometryRef().Intersects(gf):
intsec = True
ds_query = None
ds_footp = None
for fl in glob.glob(fn + '.*'):
os.remove(fl)
os.remove('B01.jp2')
clip_ok = intsec
if clip_ok:
g.write('Adding date ' + datestr[0] + ' for tile ' + tile + ' with ' + str(info['dataCoveragePercentage']) + '% data coverage and ' + str(info['cloudyPixelPercentage']) + '% cloudy pixels\n')
valid_dates.append(datestr)
else:
g.write('Rejecting date ' + datestr[0] + ' for tile ' + tile + ' with ' + str(
info['dataCoveragePercentage']) + '% data coverage and ' + str(
info['cloudyPixelPercentage']) + '% cloudy pixels (valid footprint does not intersect query)\n')
else:
g.write('Rejecting date ' + datestr[0] + ' for tile ' + tile + ' with ' + str(
info['dataCoveragePercentage']) + '% data coverage and ' + str(
info['cloudyPixelPercentage']) + '% cloudy pixels\n')
cd += datetime.timedelta(1)
if len(valid_dates) > 0:
print "Products have been found!"
for dt in valid_dates:
f.write('python getSentinel2.py ' + opt_str + ' ' + '--sequence ' + str(dt[1]) + ' ' + dt[0] + ' ' + tile + '\n')
f.close()
g.close()
DEVN.close()
if platform.system() == 'Linux':
st = os.stat('S2DownloadScript.sh')
os.chmod('S2DownloadScript.sh', st.st_mode | stat.S_IEXEC)
if __name__ == '__main__':
if len(sys.argv) < 2:
sys.exit(
'Usage: python genS2DownloadScript.py [-o <output-dir>] [-c] [-r] [-m <bufsize>] [--min-coverage <perc>] [--max-clouds <perc>] [--scan-footprint] <shapefile OR colon separated tile numbers> <start date> <end date>')
main(sys.argv[1:])