An error occurred while loading the file. Please try again.
-
Heraut Louis authored9dd131fa
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
import warnings
from TimeSeries.s2sen2cor import *
import planetary_computer as PC
from pystac_client import Client
import rasterio
import rasterio.mask
from osgeo import ogr
from pyproj import Transformer as T
from shapely.geometry import Polygon
import tqdm
def fetch(shp, dt, output_fld, band_list=None):
ds = ogr.Open(shp)
ly = ds.GetLayer(0)
xt = ly.GetExtent()
shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))
qry_srs = 4326
i_bbox = [xt[0], xt[2], xt[1], xt[3]]
if shp_srs == qry_srs:
bbox = i_bbox
else:
tr = T.from_crs(shp_srs,qry_srs,always_xy=True)
lon_min, lat_min = tr.transform(xt[0], xt[2])
lon_max, lat_max = tr.transform(xt[1], xt[3])
bbox = [lon_min, lat_min, lon_max, lat_max]
api = Client.open('https://planetarycomputer.microsoft.com/api/stac/v1', modifier=PC.sign_inplace)
res = api.search(collections="sentinel-2-l2a", bbox=bbox, datetime=dt)
lst = ['B02','B03','B04','B05','B06','B07','B08','B8A','B11','B12','SCL']
if band_list is not None:
lst = band_list
prg = tqdm.tqdm(total=len(res.item_collection()) * len(lst), desc="Fetching from Planetary")
for item in res.items():
with rasterio.open(item.assets['B02'].href) as ds:
img_srs = ds.crs.to_epsg()
if shp_srs == img_srs:
o_bbox = [20 * round(x / 20) for x in i_bbox]
else:
tr = T.from_crs(shp_srs, img_srs, always_xy=True)
x_min, y_min = tr.transform(i_bbox[0], i_bbox[1])
x_max, y_max = tr.transform(i_bbox[2], i_bbox[3])
o_bbox = [20 * round(x / 20) for x in [x_min, y_min, x_max, y_max]]
for a in lst:
ofn = os.path.join(output_fld, '/'.join(item.assets[a].get_absolute_href().split('?')[0].split('/')[10:]))
os.makedirs(os.path.dirname(ofn),exist_ok=True)
with rasterio.open(item.assets[a].href) as img:
out_img, out_geot = rasterio.mask.mask(img,
[Polygon.from_bounds(o_bbox[0], o_bbox[1],
o_bbox[2], o_bbox[3])],
crop=True)
out_meta = img.meta
out_meta.update({"driver": "GTiff",
"height": out_img.shape[1],
"width": out_img.shape[2],
"transform": out_geot})
with rasterio.open(ofn, "w", **out_meta) as dest:
dest.write(out_img)
prg.update()
prg.close()
if band_list is None:
return S2PlaneteryPipeline(output_fld)
else:
warnings.warn("Queried for a non-default band list. Skipping pipeline setup.")
return
class S2PlanetaryTilePipeline(S2Sen2CorTilePipeline):
NAME = 'S2-L2A-SEN2COR-PLANETARY'
PTRN_10m = ['GRANULE/*/IMG_DATA/R10m/*_B02_10m.tif',
'GRANULE/*/IMG_DATA/R10m/*_B03_10m.tif',
'GRANULE/*/IMG_DATA/R10m/*_B04_10m.tif',
'GRANULE/*/IMG_DATA/R10m/*_B08_10m.tif']
PTRN_20m = ['GRANULE/*/IMG_DATA/R20m/*_B05_20m.tif',
'GRANULE/*/IMG_DATA/R20m/*_B06_20m.tif',
'GRANULE/*/IMG_DATA/R20m/*_B07_20m.tif',
'GRANULE/*/IMG_DATA/R20m/*_B8A_20m.tif',
'GRANULE/*/IMG_DATA/R20m/*_B11_20m.tif',
'GRANULE/*/IMG_DATA/R20m/*_B12_20m.tif']
PTRN_msk = ['GRANULE/*/IMG_DATA/R20m/*_SCL_20m.tif']
MERG_msk = ['min']
PTRN_ful = PTRN_10m[0:3] + PTRN_20m[0:3] + [PTRN_10m[3]] + PTRN_20m[3:]
class S2PlaneteryPipeline(S2Sen2CorPipeline):
S2TilePipeline = S2PlanetaryTilePipeline
_check = S2TilePipeline._check
_tile_id = S2TilePipeline._tile_id