computeLatestClearObservation.sh 2.72 KB
Newer Older
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
#!/bin/bash
#
# Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
#
# This file is part of Let-it-snow (LIS)
#
#     https://gitlab.orfeo-toolbox.org/remote_modules/let-it-snow
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

# Script to make a composite of the latest clear-sky observations over the past nDays since now
# usage: ./computeLatestClearObservation.sh tileId inputDir outputDir nDays
# example: ./computeLatestClearObservation.sh 31TCH /work/OT/siaa/Theia/Neige/Theia/ /tmp/ 20
# output: geotiff file with colormap cyan=snow, grey=no-snow, white=cloud, black=no-data

# tile ID (ex. 31TCH)
tileId=$1
# Input directory where to look for the products
inputDir=$2
# Output directory
outputDir=$3
# Number of days to go back
nDays=$4
# Input directory containing the snow products 
#inputDir="/datalake/L2B-SNOW/${tileId}" # faster than /datalake/L2B-SNOW
#inputDir="/work/OT/siaa/Theia/Neige/Theia/"
# Today in yyyymmdd format
today=$(date +"%Y%m%d")
# name of the output composite
outputFn="${outputDir}/LCO_${tileId}_${today}.tif"
# temporary file
tmpFn="${outputDir}/tmp_${tileId}_${today}.tif"
# get date list of the n latest days
dateList=$(for i in $(seq $nDays); do date -d "$today -$i day" +"%Y%m%d"; done)
# search matching products in the tile folder and store in an array
productList=($(for d in ${dateList}; do find ${inputDir} -name "SENTINEL2[A,B]_${d}*${tileId}*SNW_R2.tif"; done))
echo $productList
# test if the product list is empty
if [ -z "$productList" ]
then
    # no products available: exit
    echo "could not find products matching these dates: "$dateList
    exit 1
else
    # load OTB
    module load otb/7.0
    # number of products to be composited
    nProduct=${#productList[@]}
    # init band math expression
    cmd=im1b1
    # build band math expression
    for i in $(seq $((nProduct-1))); do cmd=$cmd" <= 100 ? im${i}b1 : im$((i+1))b1" ; done
    # execute band math
    otbcli_BandMath -progress true -il ${productList[@]} -out "$tmpFn" uint8 -exp "$cmd"
    # apply color table to composite
    otbcli_ColorMapping -progress true -in "$tmpFn" -out "$outputFn""?&gdal:co:COMPRESS=DEFLATE" uint8 -method.custom.lut "/work/OT/siaa/Theia/hpc_scripts/LIS_SEB_style_OTB.txt"
    # delete temporary file
    rm "$tmpFn"
fi