forked from DaChro/get-sen2-py
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sen2_composites.py
72 lines (55 loc) · 3.25 KB
/
sen2_composites.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
# This script searches for the atmospherically corrected (L2A) Sentinel-2 files in a folder and creates
# 4-Band CIR composites from the 10m resolution files.
# NOTE: If Level is set to '2A' (default), this script assumes the directory structure created by sen2cor
# (which creates an additional subfolder "R10m") -> see variable bands_path. Otherwise it assumes
# the standard directory structure of Sentinel-2 data.
from sys import argv
from osgeo import gdal
import os
def create_composites(imgfolder, Level = '2A'):
# when directly run from command line, check whether composites should be created from L2A or L1C images
# when run as imported module in another script, Level can be set by that script but is not required
# create folder for composites
try:
if not os.path.exists(os.path.join(imgfolder,'CIR_composites')):
os.makedirs(os.path.join(imgfolder,'CIR_composites'))
except OSError:
print('Error: Creating directory', os.path.join(imgfolder,'CIR_composites'))
# list unzipped folders
if Level == '2A':
subfolders = [s for s in os.listdir(imgfolder) if "L2A" in s]
print("Creating 4-band CIR composites of level 2A images...")
else:
subfolders = [s for s in os.listdir(imgfolder) if "L1C" in s]
print("Creating 4-band CIR composites of level 1C images...")
# in each folder find CIR bands and create composites
for folder in subfolders:
# navigate through GRANULE folder to band files
granule_folder = os.listdir(os.path.join(imgfolder, folder, 'GRANULE'))[0]
if Level == '2A':
bands_path = os.path.join(imgfolder, folder, 'GRANULE', granule_folder, 'IMG_DATA', 'R10M')
else:
bands_path = os.path.join(imgfolder, folder, 'GRANULE', granule_folder, 'IMG_DATA')
# create list of paths to the CIR bands, search by name to be independent from number of files in folder
CIR_bandpaths = [os.path.join(bands_path, [s for s in os.listdir(bands_path) if "B02" in s][0]),
os.path.join(bands_path, [s for s in os.listdir(bands_path) if "B03" in s][0]),
os.path.join(bands_path, [s for s in os.listdir(bands_path) if "B04" in s][0]),
os.path.join(bands_path, [s for s in os.listdir(bands_path) if "B08" in s][0])]
# create virtual layer
outvrt = '/vsimem/composite.vrt' # /in-memory virtual directory
# create subfolder for current image composite
composite_path = os.path.join(imgfolder, 'CIR_composites', folder)
try:
if not os.path.exists(composite_path):
os.makedirs(composite_path)
except OSError:
print('Error: Creating directory. ' + composite_path)
# create composite in composite folder
composite_name = granule_folder.split("_")[3] + '_composite.tif'
outtif = os.path.join(composite_path, composite_name)
outds = gdal.BuildVRT(outvrt, CIR_bandpaths, separate=True)
outds = gdal.Translate(outtif, outds)
print("Composite created for " + granule_folder)
print("Finished. All composites can be found in folder " + imgfolder + '\CIR_composites.')
if __name__ == "__main__":
create_composites(argv[1],argv[2])