diff --git a/aux_data/s2_tiling_grid.dbf b/aux_data/s2_tiling_grid.dbf index 4bcbbfd9c076763cb1f7ce4a1c69180d10be7bc9..f63b669d6e0a06343628c75c77218d3d022b22a9 100644 Binary files a/aux_data/s2_tiling_grid.dbf and b/aux_data/s2_tiling_grid.dbf differ diff --git a/aux_data/s2_tiling_grid.qpj b/aux_data/s2_tiling_grid.qpj new file mode 100644 index 0000000000000000000000000000000000000000..5fbc831e743348854289f8b0cacd8b553b1262cc --- /dev/null +++ b/aux_data/s2_tiling_grid.qpj @@ -0,0 +1 @@ +GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]] diff --git a/aux_data/s2_tiling_grid.shp b/aux_data/s2_tiling_grid.shp index 8ff8bc345f929e66271e8bd203fb09fa665baaa5..ea3499483e8c7320492d5e4d360203b58a0ee4c0 100644 Binary files a/aux_data/s2_tiling_grid.shp and b/aux_data/s2_tiling_grid.shp differ diff --git a/aux_data/s2_tiling_grid.shx b/aux_data/s2_tiling_grid.shx index 57eb59e1a36511c0c82a5c34505520c7fbf8c707..5b8296e8e6cc35a1f3b585244e1579677e87e824 100644 Binary files a/aux_data/s2_tiling_grid.shx and b/aux_data/s2_tiling_grid.shx differ diff --git a/genS2DownloadScript.py b/genS2DownloadScript.py index 31e656435b0d29ce5ca8d6022ccdee04e0fd8026..2cc1175642bcf55c8ac2a1dc1075dbf80ba6bd2f 100644 --- a/genS2DownloadScript.py +++ b/genS2DownloadScript.py @@ -19,12 +19,14 @@ def findtiles(shp): clip_ly = clip.GetLayer(0) tiles = [] - ovl = [] 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() @@ -33,18 +35,22 @@ def findtiles(shp): 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) - ovl.append(intsect.GetArea() / g.GetArea()) - else: - ovl[tiles.index(tl)] += intsect.GetArea() / g.GetArea() - idx = sorted(range(len(ovl)), key=lambda k: ovl[k], reverse=True) + 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(ovl) and tot < 1.0: - tot += ovl[idx[i]] + while i < len(dovl) and tot < 1.0: + tot += dovl[idx[i]] i += 1 final_tiles = [tiles[k] for k in idx[0:i]]