Commit 18e89f9cffe67af82185f1d81f4874e7ce82779f

Authored by Germain Salgues
2 parents 9c7ec72d 17870bc4
Exists in master

Merge branch 'release/1.4'

CHANGELOG.md
1 1 # Change Log
2 2 All notable changes to LIS will be documented in this file.
3 3  
4   -## [Unreleased]
  4 +## [1.4] - 2018-02-14
  5 +
  6 +### Added
  7 +- Experimental pass1_5 function implementing the removal of snow areas inside initial cloud mask (doughnuts)
  8 +- Experimental new application to run and evaluate an annual snow map computation from a timeserie of S2 and/or L8 snow products
  9 +- Added fclear_lim parameter minimum percentage of clear pixels in an elevation band
  10 + (default value 0.1) used to compute the snow line.
  11 +- Added option to disable vector generation
  12 +- Added options to use and manage gdal_trace_outline instead of gdal_polygonize
  13 +
  14 +### Changed
  15 +- Changed default value for parameter red_darkcloud to 300 to reduce cloud sensitivity
  16 +- Changed all_cloud_mask.tif, it now include the thin clouds (in accordance with ATBD)
  17 +- Fixed method compute_percent is fix when image is empty or filled with nodata
  18 +- Fixed zs condition to trigger properly pass2 (in accordance with ATBD)
  19 +- Changed zs computation in pass2 to considers the full image imprint (including nodata pixels)
  20 +- Changed cloud mask refinement is not apply during pass1 to improve snow line accuracy during pass2
  21 +- Updated build_json.py to handle boolean parameters
  22 +- Updated build_json.py to handle new lis input parameters
  23 +
  24 +## [1.3.1] - 2017-11-23
  25 +
  26 +### Hotfix
  27 +- Fix the intermediate data format (used 1 bit instead of type uint8)
  28 +
  29 +## [1.3] - 2017-11-02
5 30  
6 31 ### Added
7 32 - Use gdal_trace_outline from the gina-alaska package instead of gdal_polygonize if available
... ...
CMakeLists.txt
... ... @@ -40,9 +40,6 @@ find_package(PythonInterp REQUIRED)
40 40 find_package( PythonLibs 2.7 REQUIRED)
41 41 include_directories( ${PYTHON_INCLUDE_DIRS} )
42 42  
43   -find_package( Boost COMPONENTS python REQUIRED )
44   -include_directories( ${Boost_INCLUDE_DIR} )
45   -
46 43 # Link to the Orfeo ToolBox
47 44 # LIS required OTB 6.0 which provides patch regarding management of 1 byte tiff image)
48 45 SET(OTB_MIN_VERSION "6.0.0")
... ...
README.md
1 1 # Let-it-snow
2 2 ## Synopsis
3 3  
4   -This code is a Python/OTB version of the snow cover extent detection algorithm for Sentinel-2 and Landsat-8 data.
  4 +This code implements the snow cover extent detection algorithm LIS (Let It Snow) for Sentinel-2, Landsat-8 and SPOT4-Take5 data.
5 5  
6 6 The algorithm documentation with examples is available here:
7 7  
... ... @@ -53,11 +53,6 @@ pixel_value & 00000101
53 53 * 100: Snow
54 54 * 205: Cloud including cloud shadow
55 55 * 254: No data
56   - * Type field:
57   - * no-snow
58   - * snow
59   - * cloud
60   - * no-data
61 56  
62 57 ## Data set example
63 58  
... ... @@ -69,23 +64,35 @@ Code to generate the snow cover extent product on Theia platform.
69 64  
70 65 ## Installation
71 66  
  67 +LIS processing chain uses CMake (http://www.cmake.org) for building from source.
  68 +
72 69 ### Dependencies
73 70  
74   -lis dependencies:
  71 +Following a summary of the required dependencies:
75 72  
76   -GDAL >=2.0
77   -OTB >= 6.0
78   -Boost-Python
79   -Python interpreter >= 2.7
80   -Python libs >= 2.7
81   -Python packages:
82   -numpy
83   -lxml
84   -matplotlib
  73 +* GDAL >=2.0
  74 +* OTB >= 6.2
  75 +* Python interpreter >= 2.7
  76 +* Python libs >= 2.7
  77 +* Python packages:
  78 +* numpy
  79 +* lxml
  80 +* matplotlib
85 81  
86 82 GDAL itself depends on a number of other libraries provided by most major operating systems and also depends on the non standard GEOS and PROJ4 libraries. GDAl- Python bindings are also required
87 83  
88   -Python package dependencies: sys, subprocess, glob, os, json, gdal
  84 +Python package dependencies:
  85 +
  86 +* sys
  87 +* subprocess
  88 +* glob
  89 +* os
  90 +* json
  91 +* gdal
  92 +
  93 +Optional dependencies:
  94 +
  95 +* gdal_trace_outline can be used alternatively to gdal_polygonize.py to generate the vector layer. It requires to install [dans-gdal-scripts utilities](https://github.com/gina-alaska/dans-gdal-scripts).
89 96  
90 97 ### Installing from the source distribution
91 98  
... ... @@ -147,7 +154,7 @@ Do not modify these folders.
147 154  
148 155 ## Contributors
149 156  
150   -Manuel Grizonnet (CNES), Simon Gascoin (CNRS/CESBIO), Tristan Klempka (CNES)
  157 +Manuel Grizonnet (CNES), Simon Gascoin (CNRS/CESBIO), Tristan Klempka (CNES), Germain Salgues (Magellium)
151 158  
152 159 ## License
153 160  
... ...
app/build_json.py
... ... @@ -11,10 +11,15 @@ conf_template = {"general":{"pout":"",
11 11 "nodata":-10000,
12 12 "ram":1024,
13 13 "nb_threads":1,
14   - "generate_vector":False,
15 14 "preprocessing":False,
16 15 "log":True,
17   - "multi":1},
  16 + "multi":1,
  17 + "target_resolution":-1},
  18 + "vector":{"generate_vector":True,
  19 + "generate_intermediate_vectors":False,
  20 + "use_gdal_trace_outline":True,
  21 + "gdal_trace_outline_dp_toler":0,
  22 + "gdal_trace_outline_min_area":0},
18 23 "inputs":{"green_band":{"path": "",
19 24 "noBand": 1},
20 25 "red_band":{"path": "",
... ... @@ -29,14 +34,20 @@ conf_template = {"general":{"pout":"",
29 34 "ndsi_pass2":0.15,
30 35 "red_pass2":40,
31 36 "fsnow_lim":0.1,
  37 + "fclear_lim":0.1,
32 38 "fsnow_total_lim":0.001},
33 39 "cloud":{"shadow_in_mask":64,
34 40 "shadow_out_mask":128,
35 41 "all_cloud_mask":1,
36 42 "high_cloud_mask":32,
37 43 "rf":12,
38   - "red_darkcloud":500,
39   - "red_backtocloud":100}}
  44 + "red_darkcloud":300,
  45 + "red_backtocloud":100,
  46 + "strict_cloud_mask":False,
  47 + "rm_snow_inside_cloud":False,
  48 + "rm_snow_inside_cloud_dilation_radius":1,
  49 + "rm_snow_inside_cloud_threshold":0.85}}
  50 +
40 51  
41 52 ### Mission Specific Parameters ###
42 53 S2_parameters = {"multi":10,
... ... @@ -88,6 +99,14 @@ mission_parameters = {"S2":S2_parameters,\
88 99 "LANDSAT8":L8_parameters,\
89 100 "Take5":Take5_parameters}
90 101  
  102 +def str2bool(v):
  103 + if v.lower() in ('yes', 'true', 't', 'y', '1'):
  104 + return True
  105 + elif v.lower() in ('no', 'false', 'f', 'n', '0'):
  106 + return False
  107 + else:
  108 + raise argparse.ArgumentTypeError('Boolean value expected.')
  109 +
91 110 def findFiles(folder, pattern):
92 111 """ Search recursively into a folder to find a patern match
93 112 """
... ... @@ -119,7 +138,7 @@ def read_product(inputPath, mission):
119 138 if result:
120 139 conf_json["inputs"]["dem"] = result[0]
121 140 else:
122   - logging.warning("No DEM found!")
  141 + logging.warning("No DEM found within product!")
123 142  
124 143 conf_json["cloud"]["shadow_in_mask"] = params["shadow_in_mask"]
125 144 conf_json["cloud"]["shadow_out_mask"] = params["shadow_out_mask"]
... ... @@ -147,14 +166,16 @@ def main():
147 166 group_general.add_argument("-nodata", type=int)
148 167 group_general.add_argument("-ram", type=int)
149 168 group_general.add_argument("-nb_threads", type=int)
150   - #group_general.add_argument("-generate_vector", type=bool)
151   - #group_general.add_argument("-preprocessing", type=bool)
152   - #group_general.add_argument("-log", type=bool)
  169 + group_general.add_argument("-generate_vector", type=str2bool, help="true/false")
  170 + group_general.add_argument("-preprocessing", type=str2bool, help="true/false")
  171 + group_general.add_argument("-log", type=str2bool, help="true/false")
153 172 group_general.add_argument("-multi", type=float)
  173 + group_general.add_argument("-target_resolution", type=float)
154 174  
155 175  
156   - group_snow = parser.add_argument_group('inputs', 'input files')
157   - group_general.add_argument("-dem", help="dem file path, to use for processing the input product")
  176 + group_inputs = parser.add_argument_group('inputs', 'input files')
  177 + group_inputs.add_argument("-dem", help="dem file path, to use for processing the input product")
  178 + group_inputs.add_argument("-cloud_mask", help="cloud mask file path")
158 179  
159 180 group_snow = parser.add_argument_group('snow', 'snow parameters')
160 181 group_snow.add_argument("-dz", type=int)
... ... @@ -173,6 +194,7 @@ def main():
173 194 group_cloud.add_argument("-rf", type=int)
174 195 group_cloud.add_argument("-red_darkcloud", type=int)
175 196 group_cloud.add_argument("-red_backtocloud", type=int)
  197 + group_cloud.add_argument("-strict_cloud_mask", type=str2bool, help="true/false")
176 198  
177 199 args = parser.parse_args()
178 200  
... ... @@ -194,21 +216,33 @@ def main():
194 216  
195 217 jsonData["general"]["pout"] = outputPath
196 218  
197   - # Overide parameters for group general
  219 + # Override parameters for group general
198 220 if args.nodata:
199 221 jsonData["general"]["nodata"] = args.nodata
  222 + if args.preprocessing is not None:
  223 + jsonData["general"]["preprocessing"] = args.preprocessing
  224 + if args.generate_vector is not None:
  225 + jsonData["vector"]["generate_vector"] = args.generate_vector
  226 + if args.log is not None:
  227 + jsonData["general"]["log"] = args.log
200 228 if args.ram:
201 229 jsonData["general"]["ram"] = args.ram
202 230 if args.nb_threads:
203 231 jsonData["general"]["nb_threads"] = args.nb_threads
204 232 if args.multi:
205 233 jsonData["general"]["multi"] = args.multi
  234 + if args.target_resolution:
  235 + jsonData["general"]["target_resolution"] = args.target_resolution
206 236  
207   - # Overide dem location
  237 + # Override dem location
208 238 if args.dem:
209 239 jsonData["inputs"]["dem"] = os.path.abspath(args.dem)
  240 + logging.warning("Using optional external DEM!")
  241 + # Override cloud mask location
  242 + if args.cloud_mask:
  243 + jsonData["inputs"]["cloud_mask"] = os.path.abspath(args.cloud_mask)
210 244  
211   - # Overide parameters for group snow
  245 + # Override parameters for group snow
212 246 if args.dz:
213 247 jsonData["snow"]["dz"] = args.dz
214 248 if args.ndsi_pass1:
... ... @@ -224,7 +258,7 @@ def main():
224 258 if args.fsnow_total_lim:
225 259 jsonData["snow"]["fsnow_total_lim"] = args.fsnow_total_lim
226 260  
227   - # Overide parameters for group cloud
  261 + # Override parameters for group cloud
228 262 if args.shadow_in_mask:
229 263 jsonData["cloud"]["shadow_in_mask"] = args.shadow_in_mask
230 264 if args.shadow_out_mask:
... ... @@ -239,6 +273,12 @@ def main():
239 273 jsonData["cloud"]["red_darkcloud"] = args.red_darkcloud
240 274 if args.red_backtocloud:
241 275 jsonData["cloud"]["red_backtocloud"] = args.red_backtocloud
  276 + if args.strict_cloud_mask:
  277 + jsonData["cloud"]["strict_cloud_mask"] = args.strict_cloud_mask
  278 +
  279 + if not jsonData["inputs"].get("dem"):
  280 + logging.error("No DEM found!")
  281 + return 1
242 282  
243 283 jsonFile = open(os.path.join(outputPath, "param_test.json"), "w")
244 284 jsonFile.write(json.dumps(jsonData, indent=4))
... ...
app/run_snow_annual_map.py 0 → 100644
... ... @@ -0,0 +1,70 @@
  1 +#!/usr/bin/env python
  2 +
  3 +import sys
  4 +import os.path as op
  5 +import json
  6 +import logging
  7 +from s2snow import snow_annual_map_evaluation
  8 +
  9 +VERSION = "0.1.0"
  10 +
  11 +
  12 +def show_help():
  13 + """Show help of the run_snow_annual_map script"""
  14 + print "This script is used to run the snow annual map " \
  15 + + "module that compute snow coverage onto a given date range"
  16 + print "Usage: python run_snow_annual_map.py param.json"
  17 + print "python run_snow_annual_map.py version to show version"
  18 + print "python run_snow_annual_map.py help to show help"
  19 +
  20 +
  21 +def show_version():
  22 + print VERSION
  23 +
  24 +# ----------------- MAIN ---------------------------------------------------
  25 +
  26 +
  27 +def main(argv):
  28 + """ main script of snow extraction procedure"""
  29 +
  30 + json_file = argv[1]
  31 +
  32 + # Load json_file from json files
  33 + with open(json_file) as json_data_file:
  34 + data = json.load(json_data_file)
  35 +
  36 + pout = data.get("path_out")
  37 + log = data.get("log", True)
  38 +
  39 + if log:
  40 + sys.stdout = open(op.join(pout, "stdout.log"), 'w')
  41 + sys.stderr = open(op.join(pout, "stderr.log"), 'w')
  42 +
  43 + # Set logging level and format.
  44 + logging.basicConfig(stream=sys.stdout, level=logging.INFO, \
  45 + format='%(asctime)s - %(filename)s:%(lineno)s - %(levelname)s - %(message)s')
  46 + logging.info("Start run_snow_annual_map.py")
  47 + logging.info("Input args = " + json_file)
  48 +
  49 + # Run the snow detector
  50 + snow_annual_map_evaluation_app = snow_annual_map_evaluation.snow_annual_map_evaluation(data)
  51 + snow_annual_map_evaluation_app.run()
  52 +
  53 + if data.get("run_l8_evaluation", False):
  54 + snow_annual_map_evaluation_app.run_evaluation()
  55 +
  56 + if data.get("run_modis_comparison", False):
  57 + snow_annual_map_evaluation_app.compare_modis()
  58 +
  59 + logging.info("End run_snow_annual_map.py")
  60 +
  61 +if __name__ == "__main__":
  62 + if len(sys.argv) != 2:
  63 + show_help()
  64 + else:
  65 + if sys.argv[1] == "version":
  66 + show_version()
  67 + elif sys.argv[1] == "help":
  68 + show_help()
  69 + else:
  70 + main(sys.argv)
... ...
app/run_snow_detector.py
... ... @@ -6,12 +6,13 @@ import json
6 6 import logging
7 7 from s2snow import snow_detector
8 8  
9   -VERSION = "1.3.1"
  9 +VERSION = "1.4"
10 10  
11 11  
12 12 def show_help():
13 13 """Show help of the run_snow_detector script"""
14   - print "This script is used to run the snow detector module that compute snow mask using OTB applications on Spot/LandSat/Sentinel-2 products from theia platform"
  14 + print "This script is used to run the snow detector module that compute snow mask" \
  15 + + " using OTB applications on Spot/LandSat/Sentinel-2 products from theia platform"
15 16 print "Usage: python run_snow_detector.py param.json"
16 17 print "python run_snow_detector.py version to show version"
17 18 print "python run_snow_detector.py help to show help"
... ... @@ -41,13 +42,15 @@ def main(argv):
41 42 sys.stderr = open(op.join(pout, "stderr.log"), 'w')
42 43  
43 44 # Set logging level and format.
44   - logging.basicConfig(stream=sys.stdout, level=logging.INFO, format='%(asctime)s - %(filename)s:%(lineno)s - %(levelname)s - %(message)s')
  45 + logging.basicConfig(stream=sys.stdout, level=logging.INFO, \
  46 + format='%(asctime)s - %(filename)s:%(lineno)s - %(levelname)s - %(message)s')
45 47 logging.info("Start run_snow_detector.py")
46 48 logging.info("Input args = " + json_file)
47 49  
48 50 # Run the snow detector
49   - sd = snow_detector.snow_detector(data)
50   - sd.detect_snow(2)
  51 + snow_detector_app = snow_detector.snow_detector(data)
  52 + snow_detector_app.detect_snow(2)
  53 + logging.info("End run_snow_detector.py")
51 54  
52 55 if __name__ == "__main__":
53 56 if len(sys.argv) != 2:
... ...
doc/actions.org
... ... @@ -1,94 +0,0 @@
1   -* CES-Neige (Suivi) [2015-11-23 lun.]
2   -** Validation S2:
3   - - premier test
4   - - récuperer image
5   -** validation landsat
6   - - sur 57 tuiles landsat
7   - - qualitatif
8   - - dégrossir
9   -** Taches
10   -*** TODO Commenter le code et ajouter des tests
11   -*** TODO avec un mnt ign, masque (données landsat), pléiades...
12   -*** TODO montrer aux utilisateurs -> pas de trou
13   -*** TODO Produits suivants
14   - - produit de synthèse: 5 ou 10 jours le temps réel en fait intéresse plus
15   - - Pour chaque date on peut interpoler avec avant/après
16   - - Intérets pour synthèse mensuel/trimestriel/
17   -*** TODO Reprise shapefile
18   - - facon de stocker le champ dans la classe -> demander à Jordi
19   - - documenter le shapefile (annoter)
20   -*** TODO script de preprocessing récupère la résolution ->
21   -*** TODO Annoncer a PS qu'on peut lancer d'images
22   -*** TODO gif animé
23   -*** TODO kml ou shapefile?
24   -*** TODO style dans le shapefile?
25   -*** TODO Renommer les fichiers de sorties
26   -*** TODO Regarder sur 2 tuiles -> comment gérer la continuité? assembler les shapefiles
27   -*** TODO Installer sur le cluster
28   -** Liens / taches
29   -Bonjour,
30   -Pour info, l'IGN vient de mettre en ligne le MNT au pas de 5 mètres sur toute la France.
31   -Nous (labos, formation) y avons accès gratuitement, comme la plupart des couches IGN.
32   -Il faut un compte pour les télécharger (j'en ai un).
33   -Je n'ai pas testé cette couche...
34   -
35   -http://professionnels.ign.fr/rgealti
36   -
37   -cdlt
38   -JFD
39   - - Les séries Landsat-8 traitées sont là:
40   -/mnt/data/home/gascoins/Landsat8/Output-CES-Neige/N2A_France-MetropoleD0005H0001/
41   - - Et la todo-list!: https://lite6.framapad.org/p/TODOs_let-it-snow
42   -
43   -* Bibliographie sur le produit neige MODIS (Gascoin) et son gap-filling
44   - 1. Le plus couramment utilisé est le produit MOD10A1 (Terra). Il y a le même avec Aqua qui
45   - s'appelle MYD10A1. Ensuite il y a un produit de synthèse au pas de temps
46   - 8-jours qui s'appelle MOD10A2 (resp. MYD10A2).
47   - 2. Voici comment est codé le raster MOD10A1 : il y a 3 bandes, une qui donne
48   - la résence/absence de neige (comme nous), une qui donne l'albédo de la
49   - neige, une qui donne la fraction de neige dans un pixel (pas pertinent pour Sentinel-2/Landsat-8):
50   - <http://nsidc.org/data/docs/daac/modis_v5/mod10a1_modis_terra_snow_daily_global_500m_grid.gd.html>
51   - Si tu descends dans la page il y a texte pédagogique avec pleins de
52   - références. Pour une doc plus technique il y a l'ATBD ici :
53   - <http://modis.gsfc.nasa.gov/data/atbd/atbd_mod10.pdf>
54   - 3. Pour le gap-filling que j'avais implémenté, tu peux regarder ce papier sur MODIS
55   - Pyrénées : Gascoin et al. (2015) <http://www.hydrol-earth-syst-sci.net/19/2337/2015/hess-19-2337-2015.html>
56   - 4. Je me suis inspiré des papiers suivants (mais il y en a d'autres !)
57   - - Parajka and Blöschl (2008) <http://www.hydro.tuwien.ac.at/fileadmin/mediapool-hydro/Publikationen/bloeschl/2008_Parajka_WRR.pdf>
58   - - Gafurov and Bárdossy (2009) <http://www.hydrol-earth-syst-sci.net/13/1361/2009/hess-13-1361-2009.pdf>
59   -
60   -* Données de validation
61   - 1. Pléiades: on a une tristéréo Pléiades sur la zone de Bassiès (Marti et al. 2016 <http://www.the-cryosphere-discuss.net/tc-2016-11/>).
62   - 2. Sur un subset de cette zone on a déjà extrait le masque de neige à partir de la panchro (mais il y a >80% de pixels neige)
63   - 3. Sinon il faut repartir des images brutes panchro et multispectrales prises le plus proche du nadir (pour la métode avec un capteur similaire voir Bühler et al. 2015 <http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6939633>)
64   -* Tâche Tristan
65   -** Test
66   -*** TODO Test n°13
67   - - Update snow fraction value
68   -** Doc
69   -*** TODO Update Readme.md
70   - - Code example
71   - - SRTM and VRT generation explanation
72   - - Files system explanation
73   - - Contributors list update
74   -** Chaine
75   -*** TODO S2Snow.py
76   - - Preprocessing integration
77   - - Cluster working directories ?
78   - - Validation tests (odd elev and snow parameters ... etc.)
79   - - Script Version
80   - - Tag 0.1
81   -*** TODO Data format
82   - - S'inspirer de la technical note pour naming. Product id => nom fichier. penser au shape file. SEB
83   - - Codage fichier 2.5.4. Notre chaine ne gère pas l'erreur. L'eau ?
84   - - Meta donnée si pass 1 != pass 2 pixel
85   - - SEB reprendre valeur du document pour notre sortie. indépendant des calculs post processing
86   - - S'assurer quon à la bonne proj wkt
87   - - Coder shapefile coder dans un champs (table attributs) pareil que SEB Type. Colonne SEB (snow, cloud, no data).
88   - - Meta data file ZS ?
89   - - QUM pass1 != pass2 ? coarse cloud => thin cloud ?
90   - - XML s'inspirer page. pytonXML 7 champs (post processing). Variable en dur
91   -*** TODO Amélioration de la chaine
92   - - Faire tourner imgs sur l'ensemble des landsat/Sentinel/Pleiades(ortho de l'image panchro avec l'OTB)
93   - - Pleaides regarder si en HR on peut générer un masqe neige regarder publi (plus tard après synthèse)
94   - - Synthèse modis atbd => algo
doc/tex/ATBD_CES-Neige.pdf
No preview for this file type
doc/tex/ATBD_CES-Neige.tex
... ... @@ -107,7 +107,25 @@ The objective of this algorithm is to generate a snow cover extent product from
107 107  
108 108 \subsection{Development}
109 109  
110   -The algorithm prototype was developed by Simon Gascoin with insights from Olivier Hagolle in June 2015. The snow detection function and a script to run this function with an example are given in appendices \ref{par:castest} and \ref{par:s2snow} as formatted documents that includes the original Matlab code, comments, and output. The LIS chain was designed to work on any high resolution multi-spectral images from satellite sensors that include at least a channel in the visible spectrum and a channel near 1.5 µm (typically referred to as mid-infrared or ``MIR''). This initial code was ported to Python 2.7 and C++ by Manuel Grizonnet in order to make it scalable to large images using Orfeo Toolbox and GDAL. The LIS code was enhanced by Tristan Klempka during his internship at CNES. LIS currently supports SPOT-4, SPOT-5, Landsat-8 and Sentinel-2 level 2A products. The LIS code, installation documentation and configuration file examples are available in the Cesbio's gitlab: \url{http://tully.ups-tlse.fr/grizonnet/let-it-snow}.
  110 +The algorithm prototype was developed by Simon Gascoin with insights from
  111 +Olivier Hagolle in June 2015. The snow detection function and a script to run
  112 +this function with an example are given in appendices \ref{par:castest} and
  113 +\ref{par:s2snow} as formatted documents that includes the original Matlab code,
  114 +comments, and output. The LIS chain was designed to work on any high resolution
  115 +multi-spectral images from satellite sensors that include at least a channel in
  116 +the visible spectrum and a channel near 1.5 µm (typically referred to as
  117 +mid-infrared or ``MIR''). This initial code was ported to Python 2.7 and C++ by
  118 +Manuel Grizonnet in order to make it scalable to large images using Orfeo
  119 +Toolbox and GDAL.
  120 +
  121 +LIS currently supports SPOT-4, SPOT-5, Landsat-8 and
  122 +Sentinel-2 level 2A products.
  123 +
  124 +The LIS code, installation documentation and configuration file examples are
  125 +available in the Cesbio's gitlab:
  126 +\url{http://tully.ups-tlse.fr/grizonnet/let-it-snow}.
  127 +
  128 +The list of all contributors is available in the LIS source in the file README.md.
111 129  
112 130 \subsection{Limitations}
113 131  
... ... @@ -166,7 +184,7 @@ The main output is a raster image (*SEB.TIF) of the snow and cloud mask. It has
166 184 \item 254: no data
167 185 \end{itemize}
168 186  
169   -The same data are made available as a polygon shapefile of the cloud and snow cover extent (*SEB\_VEC*). Two fields of information are embedded in this file:
  187 +The same data are made available as polygons (ESRI Shapefile format) of the cloud and snow cover extent (*SEB\_VEC*). Two fields of information are embedded in this file:
170 188 \begin{itemize}
171 189 \item DN:
172 190 \begin{itemize}
... ... @@ -200,9 +218,9 @@ The other output files are rather useful for the expert evaluation and troublesh
200 218  
201 219 \subsection{Pre-processing}
202 220  
203   -In the case of Sentinel-2 the red and green bands are first resampled with the cubic method to a pixel size of 20~m by 20~m to match the resolution of the MIR band.
  221 +In the case of Sentinel-2 the red and green bands are first resampled with the cubic method to a pixel size of 20~m by 20~m to match the resolution of the SWIR band.
204 222  
205   -The DEM is also resampled to the resolution of the target product (30~m or 20~m, see Sect.~\ref{par:outputs}) using the cubic spline method that is implemented in the gdal library.
  223 +The DEM is also resampled to the resolution of the target product (30~m or 20~m, see Sect.~\ref{par:outputs}) using the cubic spline method that is implemented in the GDAL library.
206 224  
207 225 \subsection{Snow detection}\label{par:snowdetec}
208 226  
... ... @@ -275,31 +293,33 @@ After passing the pass 1 and 2 snow tests, some pixels that were originally mark
275 293 \node[below right] at (leg.north west) {Legend};
276 294 \node [blockinput, below of=MUSCATE] (DEM) {DEM};
277 295 \node [decision, below of=Level 2A product] (is cloud or shadow?) {Is cloud or shadow?};
278   - \node [decision, right of=is cloud or shadow?] (darkcloud) {Is shadow or high or bright cloud? \textcolor{red}{$r_D$}};
279   - \node [blockfinal, right of=darkcloud] (cloudfinal1) {Cloud (pass 1)};
  296 + \node [blockfinal, right of=is cloud or shadow?] (cloudfinal1) {Cloud (pass 1)};
280 297 \node [decision, below of=is cloud or shadow?] (snowtest1) {Is snow? \textcolor{red}{$n_1$}, \textcolor{red}{$r_1$}};
281 298 \node [blockfinal, left of=snowtest1] (pass1) {Snow\\(pass 1)};
282 299 \node [decision, below of=snowtest1] (snowlim) {Enough snow?\\ \textcolor{red}{$f_t$}};
283 300 \node [decision, below of=snowlim] (abovezs) {Is above snowline?};
  301 + \node [decision, right of=abovezs] (darkcloud) {Is shadow or high or bright cloud? \textcolor{red}{$r_D$}};
284 302 \node [block, left of=abovezs] (zs) {Snowline elevation \textcolor{red}{$d_z$}, \textcolor{red}{$f_s$}};
285   - \node [decision, right of=abovezs] (snowtest2) {Is snow? \textcolor{red}{$n_2$}, \textcolor{red}{$r_2$}};
  303 + \node [decision, below of=darkcloud] (snowtest2) {Is snow? \textcolor{red}{$n_2$}, \textcolor{red}{$r_2$}};
286 304 \node [blockfinal, right of=snowtest2] (pass2) {Snow\\(pass 2)};
287 305 \node [decision, below of=snowtest2] (wascloud) {Was cloud?};
288 306 \node [blockfinal, right of=wascloud] (nosnow) {No snow};
289   - \node [decision, below of=wascloud] (backtocloud) {Is dark? \textcolor{red}{$r_B$}};
290   - \node [blockfinal, left of=backtocloud] (cloudfinal) {Cloud (pass 2)};
  307 + \node [decision, below of=wascloud] (backtocloud) {Is dark?
  308 + \textcolor{red}{$r_B$}};
  309 + \node [blockfinal, right of=darkcloud] (cloudfinal2) {Cloud (pass 2)};
  310 + \node [blockfinal, left of=backtocloud] (cloudfinal) {Cloud (final)};
291 311  
292 312 % Draw edges
293 313 \path [line] (Level 2A product) -- (is cloud or shadow?);
294 314 \path [line,dashed] (MUSCATE) -- (Level 2A product);
295   - \path [line] (is cloud or shadow?) -- node[near start]{yes} (darkcloud);
  315 + \path [line] (is cloud or shadow?) -- node[near start]{yes} (cloudfinal1);
296 316 \path [line] (is cloud or shadow?) -- node[near start]{no} (snowtest1);
297   - \path [line] (darkcloud) -- node[near start]{yes}(cloudfinal1);
298   - \path [line] (darkcloud) |- node[near start]{no}(snowtest1);
299 317 \path [line] (snowtest1) -- node[near start]{yes} (pass1);
300 318 \path [line] (snowtest1) -- node[near start]{yes} (snowlim);
301 319 \path [line] (snowlim) -- node[near start]{yes} (abovezs);
302   - \path [line] (abovezs) -- node[near start]{yes} (snowtest2);
  320 + \path [line] (abovezs) -- node[near start]{yes} (darkcloud);
  321 + \path [line] (darkcloud) -- node[near start]{yes}(cloudfinal2);
  322 + \path [line] (darkcloud) -- node[near start]{no}(snowtest2);
303 323 \path [line] (snowtest2) -- node[near start]{yes} (pass2);
304 324 \path [line] (snowtest2) -- node[near start]{no} (wascloud);
305 325 \path [line,dashed] (pass1) -- (zs);
... ... @@ -329,13 +349,14 @@ The table below gives the description of the main parameters of the algorithm:
329 349 Parameter & Description & Name in the configuration file & Default value\\
330 350 \hline
331 351 \textcolor{red}{$r_f$} & Resize factor to produce the down-sampled red band & \texttt{rf} & 8 for L8 (12 for S2) \\
332   -\textcolor{red}{$r_D$} & Maximum value of the down-sampled red band reflectance to define a dark cloud pixel & \texttt{rRed\_darkcloud} & 0.650 \
  352 +\textcolor{red}{$r_D$} & Maximum value of the down-sampled red band reflectance to define a dark cloud pixel & \texttt{rRed\_darkcloud} & 0.300 \
333 353 \textcolor{red}{$n_1$} & Minimum value of the NDSI for the pass 1 snow test & \texttt{ndsi\_pass1} & 0.400\\
334 354 \textcolor{red}{$n_2$} & Minimum value of the NDSI for the pass 2 snow test & \texttt{ndsi\_pass2} & 0.150\\
335 355 \textcolor{red}{$r_1$} & Minimum value of the red band reflectance the pass 1 snow test & \texttt{rRed\_pass1} & 0.200 \\
336 356 \textcolor{red}{$r_1$} & Minimum value of the red band reflectance the pass 2 snow test & \texttt{rRed\_pass2} & 0.040 \\
337   -\textcolor{red}{$d_z$} & Minimum snow fraction in an elevation band to define $z_s$ & \texttt{fsnow\_lim} & 0.100 \
  357 +\textcolor{red}{$d_z$} & Size of elevation band in the DEM used to define $z_s$ & \texttt{dz} & 0.100 \
338 358 \textcolor{red}{$f_t$} & Minimum snow fraction in an elevation band to define $z_s$ & \texttt{fsnow\_lim} & 0.100 \\
  359 +\textcolor{red}{$fc_t$} & Minimum clear pixels fraction (snow and no-snow) in an elevation band to define $z_s$ & \texttt{fclear\_lim} & 0.100 \\
339 360 \textcolor{red}{$f_s$} & Minimum snow fraction in the image to activate the pass 2 snow test & \texttt{fsnow\_total\_lim} & 0.001 \\
340 361 \textcolor{red}{$r_B$} & Minimum value of the red band reflectance to return a non-snow pixel to the cloud mask & \texttt{rRed\_backtocloud} & 0.100 \\
341 362 \hline
... ... @@ -353,12 +374,14 @@ cases, there is a parameter &#39;multi&#39; in the json configuration file which allows
353 374 to scale reflectance parameters. For instance, for products with reflectance
354 375 between 0 and 10000 you can use
355 376  
  377 +\newpage
  378 +
356 379 \subsubsection{JSON schema of configuration file}\label{par:jsonparam}
357 380  
358 381 The JSON Schema here describes the parameter file format and provide a clear, human-
359 382 and machine-readable documentation of all the algorithm parameters. JSON schema
360 383 was generated on \href{https://jsonschema.net} with the following options (with
361   -metada and relative id).
  384 +metadata and relative id).
362 385  
363 386 \inputminted[tabsize=2, fontsize=\tiny]{js}{schema.json}
364 387  
... ... @@ -388,7 +411,7 @@ The implementation of the Sentinel-2 configuration was tested on the Sentinel-2A
388 411 \centering
389 412 \includegraphics[width=\textwidth]{./images/Sentinel2_testmontage.png}
390 413 % Sentinel2_testmontage.png: 2014x811 pixel, 72dpi, 71.05x28.61 cm, bb=0 0 2014 811
391   - \caption{The LIS snow mask from the Sentinel-2A image of 06-July-2015 (Fig.~\ref{fig:S2snow}) is superposed to an aerial orthophoto taken in August 2013 and distributed by the Institut National Information Géographique Forestière.}
  414 + \caption{The LIS snow mask from the Sentinel-2A image of 06-July-2015 (Fig.~\ref{fig:S2snow}) is superposed to an aerial image taken in August 2013 and distributed by the Institut National Information Géographique Forestière.}
392 415 \label{fig:S2snowzoom}
393 416 \end{figure}
394 417  
... ...
doc/tex/schema.json
... ... @@ -27,7 +27,7 @@
27 27 "type": "integer"
28 28 },
29 29 "red_darkcloud": {
30   - "default": 500,
  30 + "default": 300,
31 31 "description": "Maximum value of the down-sampled red band reflectance to define a dark cloud pixel.",
32 32 "id": "red_darkcloud",
33 33 "title": "The Red_darkcloud schema.",
... ... @@ -53,6 +53,34 @@
53 53 "id": "shadow_out_mask",
54 54 "title": "The Shadow_out_mask schema.",
55 55 "type": "integer"
  56 + },
  57 + "strict_cloud_mask": {
  58 + "default": false,
  59 + "description": "Option that prevent any snow detection within the initial cloud mask. (experimental)",
  60 + "id": "strict_cloud_mask",
  61 + "title": "The Strict_cloud_mask schema.",
  62 + "type": "boolean"
  63 + },
  64 + "rm_snow_inside_cloud": {
  65 + "default": false,
  66 + "description": "Trigger the experimental function discarding snow area that are inside in cloud mask.",
  67 + "id": "rm_snow_inside_cloud",
  68 + "title": "The Rm_snow_inside_cloud schema.",
  69 + "type": "boolean"
  70 + },
  71 + "rm_snow_inside_cloud_dilation_radius": {
  72 + "default": 1,
  73 + "description": "Size in pixel of the dilation radius around the snow area. (experimental)",
  74 + "id": "rm_snow_inside_cloud_dilation_radius",
  75 + "title": "The Rm_snow_inside_cloud_dilation_radius schema.",
  76 + "type": "integer"
  77 + },
  78 + "rm_snow_inside_cloud_threshold": {
  79 + "default": 0.85,
  80 + "description": "Minimum fraction of cloudy pixel in the dilated area to discard the snow area. (experimental)",
  81 + "id": "rm_snow_inside_cloud_threshold",
  82 + "title": "The rm_snow_inside_cloud_threshold schema.",
  83 + "type": "float"
56 84 }
57 85 },
58 86 "type": "object"
... ... @@ -60,13 +88,6 @@
60 88 "general": {
61 89 "id": "general",
62 90 "properties": {
63   - "generate_vector": {
64   - "default": false,
65   - "description": "Generate vector masks of the detection (pass 1 and 2 and final result).",
66   - "id": "generate_vector",
67   - "title": "The Generate_vector schema.",
68   - "type": "boolean"
69   - },
70 91 "log": {
71 92 "default": true,
72 93 "description": "Log output and error to files (std***.log).",
... ... @@ -76,7 +97,7 @@
76 97 },
77 98 "multi": {
78 99 "default": 10,
79   - "description": "Scale input paramters to map reflectance interval (parameters are provided in milli-reflectance but L2A S2 Theia product are between 0 and 10000).",
  100 + "description": "Scale input parameters to map reflectance interval (parameters are provided in milli-reflectance but L2A S2 Theia product are between 0 and 10000).",
80 101 "id": "multi",
81 102 "title": "The Multi schema.",
82 103 "type": "integer"
... ... @@ -115,9 +136,56 @@
115 136 "title": "The Ram schema.",
116 137 "type": "integer"
117 138 }
  139 + "target_resolution": {
  140 + "default": -1,
  141 + "description": "Resolution of the output SNOW products in meter (automatically use the input product resolution in case target_resolution=-1)",
  142 + "id": "target_resolution",
  143 + "title": "The target_resolution schema.",
  144 + "type": "float"
  145 + }
118 146 },
119 147 "type": "object"
120 148 },
  149 + "vector": {
  150 + "id": "vector",
  151 + "properties": {
  152 + "generate_vector": {
  153 + "default": true,
  154 + "description": "Generate vector with snow and cloud masks of the final detection.",
  155 + "id": "generate_vector",
  156 + "title": "The Generate_vector schema.",
  157 + "type": "boolean"
  158 + },
  159 + "generate_intermediate_vectors": {
  160 + "default": false,
  161 + "description": "Generate vector masks of the detection (pass 1 and 2 and final result).",
  162 + "id": "generate_intermediate_vectors",
  163 + "title": "The generate_intermediate_vectors schema.",
  164 + "type": "boolean"
  165 + },
  166 + "use_gdal_trace_outline": {
  167 + "default": true,
  168 + "description": "Generate vector mask using gdal_trace_outline from gina_tools.",
  169 + "id": "use_gdal_trace_outline",
  170 + "title": "The use_gdal_trace_outline schema.",
  171 + "type": "boolean"
  172 + },
  173 + "gdal_trace_outline_min_area": {
  174 + "default": 0,
  175 + "description": "Minimum area to keep in vector mask when using gdal_trace_outline.",
  176 + "id": "gdal_trace_outline_min_area",
  177 + "title": "The gdal_trace_outline_min_area schema.",
  178 + "type": "int"
  179 + },
  180 + "gdal_trace_outline_dp_toler": {
  181 + "default": 0,
  182 + "description": "Tolered pixel approximation in vector mask when using gdal_trace_outline (0 no approximation).",
  183 + "id": "gdal_trace_outline_dp_toler",
  184 + "title": "The gdal_trace_outline_dp_toler schema.",
  185 + "type": "int"
  186 + },
  187 + "type": "object"
  188 + },
121 189 "inputs": {
122 190 "id": "inputs",
123 191 "properties": {
... ... @@ -210,6 +278,13 @@
210 278 "title": "The Fsnow_lim schema.",
211 279 "type": "number"
212 280 },
  281 + "fclear_lim": {
  282 + "default": 0.1,
  283 + "description": "Minimum clear pixel fraction in an elevation band to define zs.",
  284 + "id": "fclear_lim",
  285 + "title": "The Fclear_lim schema.",
  286 + "type": "number"
  287 + },
213 288 "fsnow_total_lim": {
214 289 "default": 0.001,
215 290 "description": "Minimum snow fraction in the image to activate the pass 2 snow test.",
... ...
hpc/makefigureTile_lis_Sentinel2_cluster_muscate.sh
... ... @@ -2,9 +2,10 @@
2 2 #PBS -N TheiaNViz
3 3 #PBS -j oe
4 4 #PBS -l select=1:ncpus=4:mem=10gb
5   -#PBS -l walltime=00:55:00
  5 +#PBS -l walltime=00:35:00
6 6 # make output figures for a better vizualisation
7   -# qsub -v tile="29SRQ" makefigureTile_lis_Sentinel2_cluster.sh
  7 +# qsub -v tile="29SRQ",input_folder="/work/OT/siaa/Theia/Neige/test_snow_in_cloud_removal" makefigureTile_lis_Sentinel2_cluster_muscate_2.sh
  8 +# useful option: qsub -W depend=afterok:<jobid> where jobid is the job id of qsub runTile..
8 9  
9 10 # IM was compiled with openMP in hal
10 11 MAGICK_THREAD_LIMIT=4 ; export MAGICK_THREAD_LIMIT
... ... @@ -15,22 +16,27 @@ export MAGICK_MAP_LIMIT MAGICK_MEMORY_LIMIT MAGICK_AREA_LIMIT
15 16 ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=4 ; export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS
16 17  
17 18 # input folder: LIS products path
18   -pin="/work/OT/siaa/Theia/Neige/output_muscate_v2/"
  19 +# Tile to process
  20 +if [ -z $input_folder ]; then
  21 + echo "input_folder is not set, exit"
  22 + exit
  23 +fi
  24 +pin=$input_folder
19 25  
20 26 # output folder: LIS figure path
21   -pout="/work/OT/siaa/Theia/Neige/output_muscate_v2/figures/"
  27 +pout=$pin"/figures/"
22 28  
23 29 # load otb
24 30 module load otb
25 31  
26 32 # Tile to process
27 33 if [ -z $tile ]; then
28   - echo "tile is not set, using 31TCH"
29   - tile=31TCH
  34 + echo "tile is not set, exit"
  35 + exit
30 36 fi
31 37  
32 38 # export colored mask in tif using otb
33   -for img in $(find $pin/T$tile/ -name *SEB.TIF)
  39 +for img in $(find $pin/T$tile/SENTINEL2A*T${tile}*D_V1-0/ -name *SEB.TIF)
34 40 do
35 41 echo $img
36 42 tiledate=$(basename $(dirname $(dirname $img)))
... ... @@ -44,16 +50,20 @@ do
44 50 echo $pout2
45 51 mkdir -p $pout2
46 52 imgout=$pout2/$labd.tif
47   - otbcli_ColorMapping -progress false -in $img -out $imgout uint8 -method.custom.lut /work/OT/siaa/Theia/hpc_scripts/LIS_SEB_style_OTB.txt
  53 + if [ -f $imgout ]; then
  54 + echo "skip $imgout (already exists)"
  55 + else
  56 + otbcli_ColorMapping -progress false -in $img -out $imgout uint8 -method.custom.lut /work/OT/siaa/Theia/hpc_scripts/LIS_SEB_style_OTB.txt
  57 + fi
48 58 # gdaldem color-relief $img /work/OT/siaa/Theia/hpc_scripts/LIS_SEB_style_v2.txt $imgout -exact_color_entry
49 59 done
50 60  
51 61 # export compo in jpg
52   -for img in $(find $pin/T$tile/ -name *COMPO.TIF)
  62 +for img in $(find $pin/T$tile/SENTINEL2A*T${tile}*D_V1-0/ -name *COMPO.TIF)
53 63 do
54 64 echo $img
55 65 tiledate=$(basename $(dirname $(dirname $img)))
56   - lab=${img:`expr index "$img" A`+1:8}
  66 + lab=${tiledate:11:8}
57 67 y=${lab:0:4}
58 68 m=${lab:4:2}
59 69 d=${lab:6:2}
... ... @@ -63,7 +73,11 @@ do
63 73 echo $pout2
64 74 mkdir -p $pout2
65 75 imgout=$pout2/$labd.jpg
66   - convert $img $imgout
  76 + if [ -f $imgout ]; then
  77 + echo "skip $imgout (already exists)"
  78 + else
  79 + convert -quiet $img $imgout
  80 + fi
67 81 done
68 82  
69 83 # make mask montage
... ...
hpc/runTile_lis_Sentinel2_cluster_muscate_anytile.sh
... ... @@ -2,17 +2,17 @@
2 2 #PBS -N TheiaNeige
3 3 #PBS -j oe
4 4 #PBS -l select=1:ncpus=1:mem=4000mb
5   -#PBS -l walltime=00:35:00
  5 +#PBS -l walltime=00:58:00
6 6 #PBS -J 1-100
7 7 #PBS -m e
8 8 # run LIS for one Sentinel-2 Level-2A tile (all available dates in the tile folder)
9 9 # specify the path to the tile folder, the path the DEM and the template configuration file (.json)
10   -# First argument is the tile name (nnccc): qsub -v tile="31TGK" runTile_lis_Sentinel2_cluster_muscate_anytile.sh
  10 +# First argument is the tile name (nnccc): qsub -v tile="31TGK",out_path="path/where/to/store/results" runTile_lis_Sentinel2_cluster_muscate_anytile.sh
11 11  
12   -# FIXME: conflict on hal between environnement variable TMPDIR used also internally by openmpi
  12 +# FIXME: conflict on hal between environnement variable TMPDIR used also internally by openmpi
13 13 # Store the directory in PBS_TMP_DIR
14 14 PBS_TMPDIR=$TMPDIR
15   -# Unset TMPDIR env variable which conflicts with openmpi
  15 +# Unset TMPDIR env variable which conflicts with openmpi
16 16 unset TMPDIR
17 17  
18 18 # Tile to process
... ... @@ -23,15 +23,21 @@ if [ -z $tile ]; then
23 23 fi
24 24 echo $tile
25 25  
  26 +if [ -z $out_path ]; then
  27 + echo "out_path is not set, exit"
  28 + exit
  29 +fi
  30 +echo $out_path
  31 +
26 32 # working directory
27 33 tmp_output_dir=$PBS_TMPDIR/TheiaNeige_Muscate_T${tile}_out/
28 34 tmp_input_dir=$PBS_TMPDIR/TheiaNeige_Muscate_T${tile}_in/
29 35  
30 36 # storage directory
31   -storage_output_dir=/work/OT/siaa/Theia/Neige/output_muscate_v2/
  37 +storage_output_dir=$out_path
32 38  
33 39 # S2 L2A products input path
34   -pin="/work/OT/muscate/prod/muscate_prod/data_production/produit/"
  40 +pin="/work/OT/siaa/Theia/S2L2A/data_production_muscate_juillet2017/L2A"
35 41  
36 42 # DEM input path
37 43 pdem="/work/OT/siaa/Theia/Neige/DEM/"
... ... @@ -82,28 +88,19 @@ mkdir -p $pout
82 88  
83 89 echo "pout" $pout
84 90  
85   -# write the config based on a template file
86   -config=$pout/$f.json
87   -cp /work/OT/siaa/Theia/hpc_scripts/lis_param_Sentinel2_template.json $config
88   -
89   -# modify only three parameters: image file, cloud file, dem file, output dir
90   -inputB11=$(find $img_input -name *FRE*B11.tif)
91   -inputB3=$(find $img_input -name *FRE*3.tif)
92   -inputB4=$(find $img_input -name *FRE*B4.tif)
93   -inputcloud=$(find $img_input -name *CLM_R2.tif)
94   -
95   -sed -i -e "s|inputB11|$inputB11|g" $config
96   -sed -i -e "s|inputB4|$inputB4|g" $config
97   -sed -i -e "s|inputB3|$inputB3|g" $config
98   -sed -i -e "s|inputcloud|$inputcloud|g" $config
99   -sed -i -e "s|inputdem|$inputdem|g" $config
100   -sed -i -e "s|outputdir|$pout|g" $config
101   -
102 91 #Load LIS modules
103 92 module load lis/develop
104 93  
105 94 #configure gdal_cachemax to speedup gdal polygonize and gdal rasterize (half of requested RAM)
106 95 export GDAL_CACHEMAX=2048
  96 +echo $GDAL_CACHEMAX
  97 +
  98 +# create config file
  99 +date ; echo "START build_json.py $config"
  100 +build_json.py -ram 2048 -dem $inputdem $img_input $pout
  101 +config=${pout}/param_test.json
  102 +echo $config
  103 +date ; echo "END build_json.py"
107 104  
108 105 # run the snow detection
109 106 date ; echo "START run_snow_detector.py $config"
... ...
hpc/runTile_lis_Sentinel2_cluster_muscate_anytile_anydate.sh
... ... @@ -5,7 +5,7 @@
5 5 #PBS -l walltime=00:55:00
6 6 # run LIS for one Sentinel-2 Level-2A tile and one date (walltime is higher)
7 7 # specify the path to the tile folder, the path the DEM and the template configuration file (.json)
8   -# First argument is the tile name (nnccc): qsub -v tile="31TCH",date="20170416" runTile_lis_Sentinel2_cluster_muscate_anytile_anydate
  8 +# First argument is the tile name (nnccc): qsub -v tile="31TCH",date="20160720",out_path="path/where/to/store/results" runTile_lis_Sentinel2_cluster_muscate_anytile_anydate
9 9 # Second argument is the date (YYYMMDD)
10 10  
11 11 # Tile to process
... ... @@ -22,15 +22,21 @@ fi
22 22  
23 23 echo $tile " on " $date
24 24  
  25 +if [ -z $out_path ]; then
  26 + echo "out_path is not set, exit"
  27 + exit
  28 +fi
  29 +echo $out_path
  30 +
25 31 # working directory
26 32 tmp_output_dir=$TMPDIR/TheiaNeige_Muscate_T${tile}_out/
27 33 tmp_input_dir=$TMPDIR/TheiaNeige_Muscate_T${tile}_in/
28 34  
29 35 # storage directory
30   -storage_output_dir=/work/OT/siaa/Theia/Neige/output_muscate_v2_debug/
  36 +storage_output_dir=$out_path
31 37  
32 38 # S2 L2A products input path
33   -pin="/work/OT/muscate/prod/muscate_prod/data_production/produit/"
  39 +pin="/work/OT/siaa/Theia/S2L2A/data_production_muscate_juillet2017/L2A"
34 40  
35 41 # DEM input path
36 42 pdem="/work/OT/siaa/Theia/Neige/DEM/"
... ... @@ -40,7 +46,7 @@ inputdem=$(find $pdem/S2__TEST_AUX_REFDE2_T${tile}_0001.DBL.DIR/ -name &quot;*ALT_R2.
40 46 echo "inputdem:" $inputdem
41 47  
42 48 # load the available product names from the tile directory
43   -array_zip=($(find $pin/SENTINEL2A_${date}*T${tile}*D* -maxdepth 2 -type f -regex ".*T${tile}.*.zip"))
  49 +array_zip=($(find $pin/SENTINEL2A_${date}*T${tile}_D* -maxdepth 2 -type f -regex ".*T${tile}.*.zip"))
44 50  
45 51 echo "array size" ${#array_zip[@]}
46 52  
... ... @@ -82,28 +88,19 @@ mkdir -p $pout
82 88  
83 89 echo "pout" $pout
84 90  
85   -# write the config based on a template file
86   -config=$pout/$f.json
87   -cp /work/OT/siaa/Theia/hpc_scripts/lis_param_Sentinel2_template.json $config
88   -
89   -# modify only three parameters: image file, cloud file, dem file, output dir
90   -inputB11=$(find $img_input -name *FRE*B11.tif)
91   -inputB3=$(find $img_input -name *FRE*3.tif)
92   -inputB4=$(find $img_input -name *FRE*B4.tif)
93   -inputcloud=$(find $img_input -name *CLM_R2.tif)
94   -
95   -sed -i -e "s|inputB11|$inputB11|g" $config
96   -sed -i -e "s|inputB4|$inputB4|g" $config
97   -sed -i -e "s|inputB3|$inputB3|g" $config
98   -sed -i -e "s|inputcloud|$inputcloud|g" $config
99   -sed -i -e "s|inputdem|$inputdem|g" $config
100   -sed -i -e "s|outputdir|$pout|g" $config
101   -
102 91 #Load LIS modules
103 92 module load lis/develop
104 93  
105 94 #configure gdal_cachemax to speedup gdal polygonize and gdal rasterize (half of requested RAM)
106 95 export GDAL_CACHEMAX=2048
  96 +echo $GDAL_CACHEMAX
  97 +
  98 +# create config file
  99 +date ; echo "START build_json.py $config"
  100 +build_json.py -ram 2048 -dem $inputdem $img_input $pout
  101 +config=${pout}/param_test.json
  102 +echo $config
  103 +date ; echo "END build_json.py"
107 104  
108 105 # run the snow detection
109 106 date ; echo "START run_snow_detector.py $config"
... ...
python/s2snow/app_wrappers.py
... ... @@ -164,6 +164,114 @@ def band_mathX(il, out, exp, ram=None, out_type=None):
164 164 else:
165 165 logging.error("Parameters il, out and exp are required")
166 166  
  167 +def compute_snow_line(img_dem, img_snow, img_cloud, dz, fsnowlim, fclearlim, \
  168 + reverse, offset, centeroffset, outhist, ram=None):
  169 + """ Create and configure the ComputeSnowLine application
  170 + using otb.Registry.CreateApplication("ComputeSnowLine")
  171 +
  172 + Keyword arguments:
  173 + TODO
  174 + """
  175 + if img_dem and img_snow and img_cloud:
  176 + logging.info("Processing ComputeSnowLine with args:")
  177 + logging.info(img_dem)
  178 + logging.info(img_snow)
  179 + logging.info(img_cloud)
  180 + logging.info(dz)
  181 + logging.info(fsnowlim)
  182 + logging.info(outhist)
  183 +
  184 + snowLineApp = otb.Registry.CreateApplication("ComputeSnowLine")
  185 + snowLineApp.SetParameterString("dem", img_dem)