Commit b40b65b1850c3ee827d7d03b8f1e85401bdd0460

Authored by Grizonnet Manuel (CNES)
2 parents e154be98 038f3910
Exists in develop and in 1 other branch master

Merge branch 'release/1.2'

.gitignore
... ... @@ -11,4 +11,5 @@ doc/tex/*aux
11 11 doc/tex/*log
12 12 doc/tex/*out
13 13 doc/tex/*toc
  14 +doc/tex/_minted*
14 15  
... ...
CHANGELOG.md 0 → 100644
... ... @@ -0,0 +1,22 @@
  1 +# Change Log
  2 +All notable changes to LIS will be documented in this file.
  3 +
  4 +## [1.2] - 2017-06-04
  5 +- add json schema to ATBD to document all parameters
  6 +- Add version of lis in atbd
  7 +- Document how to build the documentation in doc/tex directory
  8 +- Compact histogram files and copy it in LIS_PRODUCTS
  9 +- Apply autopep8 to all Python scripts to imprve code quality
  10 +- Add a changelog
  11 +## [1.1.1] - 2016-11-28
  12 +- minor update in build scripts
  13 +- change ctest launcher location
  14 +## [1.1.0] - 2016-11-28
  15 +- Change license from GPL to A-GPL
  16 +- Improvments in cmake configuration
  17 +- launch tests in separate directories
  18 +## [1.0.0] - 2016-07-06
  19 +- First released version of LIS with support with last MUSCATE format
  20 +- Support for image splitted in multiple files
  21 +- Use high clouds mask in cloud refine
  22 +- Prototype for multi-T synthesis
... ...
README.md
... ... @@ -11,32 +11,27 @@ To read more about the "Centre d'Expertise Scientifique surface enneigée" (in F
11 11  
12 12 * [Bulletin THEIA](https://www.theia-land.fr/sites/default/files/imce/BulletinTHEIA3.pdf#page=10)
13 13  
14   -The input files are Sentinel-2 or Landsat-8 level-2A products from the [Theai Land Data Centre](https://theia.cnes.fr/) or [SPOT-4/5 Take 5 level-2A products](https://spot-take5.org) and the SRTM digital elevation model.
  14 +The input files are Sentinel-2 or Landsat-8 level-2A products from the [Theai Land Data Centre](https://theia.cnes.fr/) or [SPOT-4/5 Take 5 level-2A products](https://spot-take5.org) and the SRTM digital elevation model reprojected at the same resolution as the input image.
15 15  
16   -## Code Example
  16 +##Usage
17 17  
18   -To build DEM data download the SRTM files corresponding to the study area and build the .vrt using gdalbuildvrt. Edit config.json file to activate preprocessing : Set "preprocessing" to true and set the vrt path.
  18 +Run the python script run_snow_detector.py with a json configuration file as unique argument:
19 19  
20   -The snow detection is performed in the Python script app/run_snow_detector.py.
21   -
22   -Configure PYTHONPATH environment
23 20 ```bash
24   -export PYTHONPATH=${lis-build-dir}/app/:$PYTHONPATH
  21 +python run_snow_detector.py param.json
25 22 ```
26   -Run the main python script:
  23 +The snow detection is performed in the Python script [run_snow_detector.py](app/run_snow_detector.py).
27 24  
28   -```bash
29   -run_snow_detector.py param.json
30   -```
  25 +All the parameters of the algorithm, paths to input and output data are stored in the json file. See the provided example [param_test_s2_template.json](tes/param_test_s2_template.json) file for an example.
  26 +
  27 +Moreover The JSON schema is available in the [Algorithm theoritical basis documentation](doc/tex/ATBD_CES-Neige.tex) and gives more information about the roles of these parameters.
  28 +
  29 +NB: To build DEM data download the SRTM files corresponding to the study area and build the .vrt using gdalbuildvrt. Edit config.json file to activate preprocessing : Set "preprocessing" to true and set the vrt path.
31 30  
32   -There is a Bash script in app directory which allows to set the env variable and run the script:
33 31  
34   -```bash
35   -runLis.sh param.json
36   -```
37 32 ## Products format
38 33  
39   -* COMPO: RGB composition with snow mask
  34 +* COMPO: Raster image showing the outlines of the cloud (including cloud shadow) and snow masks drawn on the RGB composition of the L2A image (bands SWIR/Red/Green).
40 35 * SNOW_ALL: Binary mask of snow and clouds.
41 36 * 1st bit: Snow mask after pass1
42 37 * 2nd bit: Snow mask after pass2
... ...
app/build_json.py 0 → 100644
... ... @@ -0,0 +1,71 @@
  1 +#!/usr/bin/env python
  2 +
  3 +import sys
  4 +import os.path as op
  5 +import json
  6 +import argparse
  7 +
  8 +def show_help():
  9 + """Show help of the build_json script for theia N2A products"""
  10 + print "This script is used to build json configuration file use then to compute snow mask using OTB applications on Spot/LandSat/Sentinel-2 products from theia platform"
  11 + print "Usage: python build_theia_json -s [landsat|s2|take5] -d image_directory -e srtm_tile -o file.json"
  12 + print "python run_snow_detector.py help to show help"
  13 +
  14 +
  15 +#----------------- MAIN ---------------------------------------------------
  16 +
  17 +def main():
  18 + """ Script to build json from theia N2A product"""
  19 +
  20 + parser = argparse.ArgumentParser(description='Build json from THEIA product')
  21 +
  22 + parser.add_argument("-s", help="select input sensors")
  23 + parser.add_argument("-d", help="input dir")
  24 + parser.add_argument("-o", help="input dir")
  25 + parser.add_argument("-do", help="input dir")
  26 +
  27 + args = parser.parse_args()
  28 +
  29 + #print(args.accumulate(args.integers))
  30 +
  31 + #Parse sensor
  32 + if (args.s == 's2'):
  33 + multi=10
  34 + #Build json file
  35 + data = {}
  36 +
  37 + data["general"]={
  38 + "pout":args.do,
  39 + "nodata":-10000,
  40 + "ram":1024,
  41 + "nb_threads":1,
  42 + "generate_vector":"false",
  43 + "preprocessing":"false",
  44 + "log":"true",
  45 + "multi":10
  46 + }
  47 +
  48 + data["cloud"]={
  49 + "shadow_mask":32,
  50 + "all_cloud_mask":1,
  51 + "high_cloud_mask":128,
  52 + "rf":12,
  53 + "red_darkcloud":500,
  54 + "red_backtocloud":100
  55 + }
  56 + data["snow"]={
  57 + "dz":100,
  58 + "ndsi_pass1":0.4,
  59 + "red_pass1":200,
  60 + "ndsi_pass2":0.15,
  61 + "red_pass2":120,
  62 + "fsnow_lim":0.1,
  63 + "fsnow_total_lim":0.001
  64 + }
  65 +
  66 + fp = open(args.o, 'w')
  67 + fp.write(json.dumps(data,indent=4, sort_keys=True))
  68 + fp.close()
  69 +
  70 +if __name__ == "__main__":
  71 + main()
... ...
app/run_snow_detector.py
... ... @@ -5,7 +5,8 @@ import os.path as op
5 5 import json
6 6 from s2snow import snow_detector
7 7  
8   -VERSION="1.1.1"
  8 +VERSION = "1.2"
  9 +
9 10  
10 11 def show_help():
11 12 """Show help of the run_snow_detector script"""
... ... @@ -14,34 +15,35 @@ def show_help():
14 15 print "python run_snow_detector.py version to show version"
15 16 print "python run_snow_detector.py help to show help"
16 17  
  18 +
17 19 def show_version():
18 20 print VERSION
19 21  
20   -#----------------- MAIN ---------------------------------------------------
  22 +# ----------------- MAIN ---------------------------------------------------
  23 +
21 24  
22 25 def main(argv):
23 26 """ main script of snow extraction procedure"""
24 27  
25   - json_file=argv[1]
  28 + json_file = argv[1]
26 29  
27   - #load json_file from json files
  30 + # Load json_file from json files
28 31 with open(json_file) as json_data_file:
29   - data = json.load(json_data_file)
30   -
  32 + data = json.load(json_data_file)
  33 +
31 34 general = data["general"]
32 35 pout = general.get("pout")
33   -
34 36 log = general.get("log", True)
  37 +
35 38 if log:
36 39 sys.stdout = open(op.join(pout, "stdout.log"), 'w')
37 40 sys.stderr = open(op.join(pout, "stderr.log"), 'w')
38   -
  41 +
39 42 sd = snow_detector.snow_detector(data)
40   -
41 43 sd.detect_snow(2)
42   -
  44 +
43 45 if __name__ == "__main__":
44   - if len(sys.argv) != 2 :
  46 + if len(sys.argv) != 2:
45 47 show_help()
46 48 else:
47 49 if sys.argv[1] == "version":
... ...
config/param_test_s2.json
... ... @@ -29,7 +29,7 @@
29 29 "shadow_mask":64,
30 30 "all_cloud_mask":0,
31 31 "high_cloud_mask":32,
32   - "rf":8,
  32 + "rf":12,
33 33 "red_darkcloud":650,
34 34 "red_backtocloud":100
35 35 },
... ...
doc/tex/ATBD_CES-Neige.pdf
No preview for this file type
doc/tex/ATBD_CES-Neige.tex
... ... @@ -30,7 +30,8 @@ pdftitle=Algorithm theoretical basis documentation for an operational snow cover
30 30 \usepackage[]{algorithm2e}
31 31 \usepackage{fancyhdr}
32 32 \usepackage{tabularx}
33   -% \usepackage{listings}
  33 +%\usepackage{listings}
  34 +\usepackage{minted}
34 35 % \lstset{language=Matlab}
35 36  
36 37 \usepackage{tikz}
... ... @@ -227,7 +228,8 @@ The snow detection (Sect.~\ref{par:snowdetec}) is performed a first time using t
227 228 The L2A cloud mask is conservative because it is computed at a coarser resolution and also because it is developed for a large range of applications. However, the detection of the snow cover is robust to a thin, transparent, cloud cover. More importantly, the L2A cloud mask tends to falsely classify the edges of the snow cover as cloud. Hence, it is possible to recover many pixels from the L2A cloud mask and reclassify them as snow or no-snow. This step is important because it substantially increases the number of observations. A pixel from the L2A cloud mask cannot be reclassified as snow or no-snow if:
228 229  
229 230 \begin{itemize}
230   - \item it is coded as ``cloud shadow'' in L2A cloud mask;
  231 + \item it is coded as ``cloud shadow'' in L2A cloud mask. Note that it can be
  232 + cloud shadows matched with a cloud or cloud shadows in the zone where clouds could be outside the image ;
231 233 \item or: it is coded as ``high altitude cloud'' (or ``cirrus'') in the L2A cloud mask;
232 234 \item or: it is not a ``dark'' cloud (see below).
233 235 \end{itemize}
... ... @@ -256,7 +258,7 @@ After passing the pass 1 and 2 snow tests, some pixels that were originally mark
256 258 minimum height=2em]
257 259 \tikzstyle{bigbox}=[inner sep=20pt]
258 260  
259   -\begin{figure}
  261 +\begin{figure}[H]
260 262  
261 263 \begin{tikzpicture}[node distance = 3.5cm, auto]
262 264 % Place nodes
... ... @@ -314,9 +316,13 @@ After passing the pass 1 and 2 snow tests, some pixels that were originally mark
314 316 \caption{Flowchart of the snow detection algorithm}
315 317 \end{figure}\label{fig:flowchart}
316 318  
317   -\subsection{Parameter values}\label{par:param}
  319 +\subsection{Parameters description}\label{par:param}
318 320  
319   -\begin{table}
  321 +\subsubsection{Main algorithm parameters}\label{par:sciparam}
  322 +
  323 +The table below gives the description of the main parameters of the algorithm:
  324 +
  325 +\begin{table}[!htbp]
320 326 \begin{center}
321 327 \begin{tabularx}{\textwidth}{|l X l l|}
322 328 \hline
... ... @@ -338,7 +344,23 @@ Parameter & Description & Name in the configuration file & Default value\\
338 344 \caption{LIS algorithm parameters description and default values.}
339 345 \end{table}\label{tab:param}
340 346  
  347 +Above default values related to reflectance are given as float values between 0
  348 +and 1. Threshold related to reflectance values follow the convention of
  349 +considering milli-reflectance as input (values between 0 and 1000) in the json
  350 +file. Some products can encode reflectance with other convention (floating
  351 +values between 0 and 1 or reflectance between 0 and 10000), to handle those
  352 +cases, there is a parameter 'multi' in the json configuration file which allows
  353 +to scale reflectance parameters. For instance, for products with reflectance
  354 +between 0 and 10000 you can use
  355 +
  356 +\subsubsection{JSON schema of configuration file}\label{par:jsonparam}
  357 +
  358 +The JSON Schema here describes the parameter file format and provide a clear, human-
  359 +and machine-readable documentation of all the algorithm parameters. JSON schema
  360 +was generated on \href{https://jsonschema.net} with the following options (with
  361 +metada and relative id).
341 362  
  363 +\inputminted[tabsize=2, fontsize=\tiny]{js}{schema.json}
342 364  
343 365 \section{Validation}\label{par:validation}
344 366  
... ...
doc/tex/README.md 0 → 100644
... ... @@ -0,0 +1,9 @@
  1 +# Let-it-snow documentation
  2 +
  3 +## Build instructions
  4 +
  5 +```bash
  6 +texi2pdf --shell-escape ATBD_CES-Neige.tex
  7 +```
  8 +
  9 +The "--shell-escape" option is required by the minted package
... ...
doc/tex/page_titre.tex
... ... @@ -14,7 +14,7 @@
14 14 { \huge \bfseries Algorithm theoritical basis documentation for an operational snow cover extent product from Sentinel-2 and Landsat-8 data (Let-it-snow)\\}
15 15 \rule{\linewidth}{0.5mm}
16 16 { \large \bfseries Simon Gascoin (CNRS/CESBIO), Manuel Grizonnet (CNES/CT/DSI), Tristan Klempka (CNES/CT/DSI)\\ }
17   -{ \large \bfseries V1.0 - \today \ }
  17 +{ \large \bfseries V1.0 (Updated for LIS 1.2) - \today \ }
18 18  
19 19 % \vspace{3cm}
20 20 % \includegraphics[width=1\textwidth]{./Images/Theia_en.png}
... ...
doc/tex/schema.json 0 → 100644
... ... @@ -0,0 +1,253 @@
  1 +{
  2 + "$schema": "http://json-schema.org/draft-04/schema#",
  3 + "id": "http://tully.ups-tlse.fr/grizonnet/let-it-snow/blob/master/test/param_test_s2_template.json",
  4 + "properties": {
  5 + "cloud": {
  6 + "id": "cloud",
  7 + "properties": {
  8 + "all_cloud_mask": {
  9 + "default": 1,
  10 + "description": "Threshold apply to Theia cloud mask to retrieve a strict cloud mask (greater than all_cloud_mask).",
  11 + "id": "all_cloud_mask",
  12 + "title": "The All_cloud_mask schema.",
  13 + "type": "integer"
  14 + },
  15 + "high_cloud_mask": {
  16 + "default": 128,
  17 + "description": "bitmask apply to retrieve high clouds for input mask",
  18 + "id": "high_cloud_mask",
  19 + "title": "The High_cloud_mask schema.",
  20 + "type": "integer"
  21 + },
  22 + "red_backtocloud": {
  23 + "default": 100,
  24 + "description": "Minimum value of the red band reflectance to return a non-snow pixel to the cloud mask.",
  25 + "id": "red_backtocloud",
  26 + "title": "The Red_backtocloud schema.",
  27 + "type": "integer"
  28 + },
  29 + "red_darkcloud": {
  30 + "default": 500,
  31 + "description": "Maximum value of the down-sampled red band reflectance to define a dark cloud pixel.",
  32 + "id": "red_darkcloud",
  33 + "title": "The Red_darkcloud schema.",
  34 + "type": "integer"
  35 + },
  36 + "rf": {
  37 + "default": 12,
  38 + "description": "Resize factor to produce the down-sampled red band (use for cloud refinement).",
  39 + "id": "rf",
  40 + "title": "The Rf schema.",
  41 + "type": "integer"
  42 + },
  43 + "shadow_in_mask": {
  44 + "default": 32,
  45 + "description": "bitmask apply to retrieve cloud shadows (for cloud inside the image).",
  46 + "id": "shadow_in_mask",
  47 + "title": "The Shadow_in_mask schema.",
  48 + "type": "integer"
  49 + },
  50 + "shadow_out_mask": {
  51 + "default": 64,
  52 + "description": "bitmask apply to retrieve cloud shadows (for cloud outside the image).",
  53 + "id": "shadow_out_mask",
  54 + "title": "The Shadow_out_mask schema.",
  55 + "type": "integer"
  56 + }
  57 + },
  58 + "type": "object"
  59 + },
  60 + "general": {
  61 + "id": "general",
  62 + "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 + "log": {
  71 + "default": true,
  72 + "description": "Log output and error to files (std***.log).",
  73 + "id": "log",
  74 + "title": "The Log schema.",
  75 + "type": "boolean"
  76 + },
  77 + "multi": {
  78 + "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).",
  80 + "id": "multi",
  81 + "title": "The Multi schema.",
  82 + "type": "integer"
  83 + },
  84 + "nb_threads": {
  85 + "default": 1,
  86 + "description": "Maximum number of threads use by the program.",
  87 + "id": "nb_threads",
  88 + "title": "The Nb_threads schema.",
  89 + "type": "integer"
  90 + },
  91 + "nodata": {
  92 + "default": -10000,
  93 + "description": "No-data value in the input L2A product.",
  94 + "id": "nodata",
  95 + "title": "The Nodata schema.",
  96 + "type": "integer"
  97 + },
  98 + "pout": {
  99 + "description": "Path to output directory.",
  100 + "id": "pout",
  101 + "title": "The Pout schema.",
  102 + "type": "string"
  103 + },
  104 + "preprocessing": {
  105 + "default": false,
  106 + "description": "Activate the extraction and resampling of the DEM.",
  107 + "id": "preprocessing",
  108 + "title": "The Preprocessing schema.",
  109 + "type": "boolean"
  110 + },
  111 + "ram": {
  112 + "default": 1024,
  113 + "description": "Maximum number of RAM memory used by the program.",
  114 + "id": "ram",
  115 + "title": "The Ram schema.",
  116 + "type": "integer"
  117 + }
  118 + },
  119 + "type": "object"
  120 + },
  121 + "inputs": {
  122 + "id": "inputs",
  123 + "properties": {
  124 + "cloud_mask": {
  125 + "description": "Input mask image (in MASK directory for Theia product).",
  126 + "id": "cloud_mask",
  127 + "title": "The Cloud_mask schema.",
  128 + "type": "string"
  129 + },
  130 + "dem": {
  131 + "description": "Input DEM with the same resolution and extent as the input image if preprocessing is deactivated.",
  132 + "id": "dem",
  133 + "title": "The Dem schema.",
  134 + "type": "string"
  135 + },
  136 + "green_band": {
  137 + "id": "green_band",
  138 + "properties": {
  139 + "noBand": {
  140 + "default": 1,
  141 + "description": "Green band number.",
  142 + "id": "noBand",
  143 + "title": "The Noband schema.",
  144 + "type": "integer"
  145 + },
  146 + "path": {
  147 + "description": "Path to input L2A image or L2A green band.",
  148 + "id": "path",
  149 + "title": "The Path schema.",
  150 + "type": "string"
  151 + }
  152 + },
  153 + "type": "object"
  154 + },
  155 + "red_band": {
  156 + "id": "red_band",
  157 + "properties": {
  158 + "noBand": {
  159 + "default": 1,
  160 + "description": "Red band number.",
  161 + "id": "noBand",
  162 + "title": "The Noband schema.",
  163 + "type": "integer"
  164 + },
  165 + "path": {
  166 + "description": "Path to input L2A image or L2A red band.",
  167 + "id": "path",
  168 + "title": "The Path schema.",
  169 + "type": "string"
  170 + }
  171 + },
  172 + "type": "object"
  173 + },
  174 + "swir_band": {
  175 + "id": "swir_band",
  176 + "properties": {
  177 + "noBand": {
  178 + "default": 1,
  179 + "description": "SWIR band number.",
  180 + "id": "noBand",
  181 + "title": "The Noband schema.",
  182 + "type": "integer"
  183 + },
  184 + "path": {
  185 + "description": "Path to input L2A image or L2A swir band.",
  186 + "id": "path",
  187 + "title": "The Path schema.",
  188 + "type": "string"
  189 + }
  190 + },
  191 + "type": "object"
  192 + }
  193 + },
  194 + "type": "object"
  195 + },
  196 + "snow": {
  197 + "id": "snow",
  198 + "properties": {
  199 + "dz": {
  200 + "default": 100,
  201 + "description": "Minimum snow fraction in an elevation band to define zs.",
  202 + "id": "dz",
  203 + "title": "The Dz schema.",
  204 + "type": "integer"
  205 + },
  206 + "fsnow_lim": {
  207 + "default": 0.1,
  208 + "description": "Minimum snow fraction in an elevation band to define zs.",
  209 + "id": "fsnow_lim",
  210 + "title": "The Fsnow_lim schema.",
  211 + "type": "number"
  212 + },
  213 + "fsnow_total_lim": {
  214 + "default": 0.001,
  215 + "description": "Minimum snow fraction in the image to activate the pass 2 snow test.",
  216 + "id": "fsnow_total_lim",
  217 + "title": "The Fsnow_total_lim schema.",
  218 + "type": "number"
  219 + },
  220 + "ndsi_pass1": {
  221 + "default": 0.4,
  222 + "description": "Minimum value of the NDSI for the pass 1 snow test.",
  223 + "id": "ndsi_pass1",
  224 + "title": "The Ndsi_pass1 schema.",
  225 + "type": "number"
  226 + },
  227 + "ndsi_pass2": {
  228 + "default": 0.15,
  229 + "description": "Minimum value of the NDSI for the pass 2 snow test.",
  230 + "id": "ndsi_pass2",
  231 + "title": "The Ndsi_pass2 schema.",
  232 + "type": "number"
  233 + },
  234 + "red_pass1": {
  235 + "default": 200,
  236 + "description": "Minimum value of the red band reflectance the pass 1 snow test.",
  237 + "id": "red_pass1",
  238 + "title": "The Red_pass1 schema.",
  239 + "type": "integer"
  240 + },
  241 + "red_pass2": {
  242 + "default": 120,
  243 + "description": "Minimum value of the red band reflectance the pass 2 snow test.",
  244 + "id": "red_pass2",
  245 + "title": "The Red_pass2 schema.",
  246 + "type": "integer"
  247 + }
  248 + },
  249 + "type": "object"
  250 + }
  251 + },
  252 + "type": "object"
  253 +}
... ...
python/s2snow/cloud_builder.py
1   -import sys, os, numpy, random
  1 +import sys
  2 +import os
  3 +import numpy
  4 +import random
2 5 import os.path as op
3 6  
4   -import gdal, gdalconst
  7 +import gdal
  8 +import gdalconst
5 9 from subprocess import call
6 10  
  11 +
7 12 def show_help():
8   - print "This script is used to create clouds on data"
9   - print "Usage: cloud_builder.py mode plaincloudthreshold randomcloudthreshold inputpath outputplaincloudpath ouputrandomcloudpath"
10   - print "Mode : 0 %plain cloud image, 1 %random cloud image, 2 both"
  13 + print "This script is used to create clouds on data"
  14 + print "Usage: cloud_builder.py mode plaincloudthreshold randomcloudthreshold inputpath outputplaincloudpath ouputrandomcloudpath"
  15 + print "Mode : 0 %plain cloud image, 1 %random cloud image, 2 both"
  16 +
  17 +
11 18 def main(argv):
12 19 mode = argv[1]
13 20 mode = int(mode)
... ... @@ -17,35 +24,38 @@ def main(argv):
17 24 input_path = argv[4]
18 25 output_plain_cloud_path = argv[5]
19 26 output_random_cloud_path = argv[6]
20   -
  27 +
21 28 dataset = gdal.Open(input_path, gdalconst.GA_ReadOnly)
22 29 wide = dataset.RasterXSize
23 30 high = dataset.RasterYSize
24   -
  31 +
25 32 if(mode == 0 or mode == 2):
26   - #build half cloud image
27   - str_exp = "idxX>"+str(int(wide*plain_cloud_threshold))+"?im1b1:205"
28   - call(["otbcli_BandMathX", "-il", input_path, "-out", output_plain_cloud_path, "-exp", str_exp])
  33 + # build half cloud image
  34 + str_exp = "idxX>" + \
  35 + str(int(wide * plain_cloud_threshold)) + "?im1b1:205"
  36 + call(["otbcli_BandMathX", "-il", input_path, "-out",
  37 + output_plain_cloud_path, "-exp", str_exp])
29 38 if(mode == 1 or mode == 2):
30   - #build random cloud image
  39 + # build random cloud image
31 40 band = dataset.GetRasterBand(1)
32 41 array = band.ReadAsArray(0, 0, wide, high)
33 42 for i in range(0, wide):
34 43 for j in range(0, high):
35 44 if(random.randint(0, 100) < random_cloud_threshold):
36   - array[i,j] = 205
37   -
38   - output_random_cloud_raster = gdal.GetDriverByName('GTiff').Create(output_random_cloud_path, wide, high, 1 ,gdal.GDT_Byte)
39   - output_random_cloud_raster.GetRasterBand(1).WriteArray(array)
  45 + array[i, j] = 205
  46 +
  47 + output_random_cloud_raster = gdal.GetDriverByName('GTiff').Create(
  48 + output_random_cloud_path, wide, high, 1, gdal.GDT_Byte)
  49 + output_random_cloud_raster.GetRasterBand(1).WriteArray(array)
40 50 output_random_cloud_raster.FlushCache()
41 51  
42 52 # georeference the image and set the projection
43 53 output_random_cloud_raster.SetGeoTransform(dataset.GetGeoTransform())
44   - output_random_cloud_raster.SetProjection(dataset.GetProjection())
45   -
  54 + output_random_cloud_raster.SetProjection(dataset.GetProjection())
  55 +
  56 +
46 57 if __name__ == "__main__":
47 58 if len(sys.argv) != 7:
48 59 show_help()
49 60 else:
50 61 main(sys.argv)
51   -
... ...
python/s2snow/cloud_removal.py
1   -import os, sys, json, multiprocessing, csv
  1 +import os
  2 +import sys
  3 +import json
  4 +import multiprocessing
  5 +import csv
2 6 import numpy as np
3 7 import matplotlib.pyplot as plot
4 8 from subprocess import call
... ... @@ -6,296 +10,468 @@ import os.path as op
6 10 import gdal
7 11 import gdalconst
8 12  
  13 +
9 14 def show_help():
10   - print "This script is used to remove clouds from snow data"
11   - print "Usage: cloud_removal.py config.json"
  15 + print "This script is used to remove clouds from snow data"
  16 + print "Usage: cloud_removal.py config.json"
  17 +
12 18  
13 19 def get_raster_as_array(raster_file_name):
14   - dataset = gdal.Open(raster_file_name, gdalconst.GA_ReadOnly)
15   - wide = dataset.RasterXSize
16   - high = dataset.RasterYSize
17   - band = dataset.GetRasterBand(1)
18   - array = band.ReadAsArray(0, 0, wide, high)
19   - return array, dataset
  20 + dataset = gdal.Open(raster_file_name, gdalconst.GA_ReadOnly)
  21 + wide = dataset.RasterXSize
  22 + high = dataset.RasterYSize
  23 + band = dataset.GetRasterBand(1)
  24 + array = band.ReadAsArray(0, 0, wide, high)
  25 + return array, dataset
  26 +
20 27  
21 28 def set_array_as_raster(array, dataset, output_path):
22   - high, wide = array.shape
23   - output = gdal.GetDriverByName('GTiff').Create(output_path, wide, high, 1 ,gdal.GDT_Byte)
24   - output.GetRasterBand(1).WriteArray(array)
25   -
26   - # georeference the image and set the projection
27   - output.SetGeoTransform(dataset.GetGeoTransform())
28   - output.SetProjection(dataset.GetProjection())
29   - output = None
  29 + high, wide = array.shape
  30 + output = gdal.GetDriverByName('GTiff').Create(
  31 + output_path, wide, high, 1, gdal.GDT_Byte)
  32 + output.GetRasterBand(1).WriteArray(array)
  33 +
  34 + # georeference the image and set the projection
  35 + output.SetGeoTransform(dataset.GetGeoTransform())
  36 + output.SetProjection(dataset.GetProjection())
  37 + output = None
  38 +
  39 +
30 40 def compute_HSmin(image_path, dem_path):
31   - array_image, dataset_image = get_raster_as_array(image_path)
32   - array_dem, dataset_dem = get_raster_as_array(dem_path)
  41 + array_image, dataset_image = get_raster_as_array(image_path)
  42 + array_dem, dataset_dem = get_raster_as_array(dem_path)
  43 +
  44 + return array_dem[array_image == 100].min()
33 45  
34   - return array_dem[array_image == 100].min()
35 46  
36 47 def compute_HSmax(image_path, dem_path):
37   - array_image, dataset_image = get_raster_as_array(image_path)
38   - array_dem, dataset_dem = get_raster_as_array(dem_path)
  48 + array_image, dataset_image = get_raster_as_array(image_path)
  49 + array_dem, dataset_dem = get_raster_as_array(dem_path)
  50 +
  51 + elev_max_nosnow = array_dem[array_image == 0].max()
  52 + return array_dem[(array_image == 100) & (
  53 + array_dem > elev_max_nosnow)].min()
39 54  
40   - elev_max_nosnow = array_dem[array_image == 0].max()
41   - return array_dem[(array_image == 100) & (array_dem > elev_max_nosnow)].min()
42 55  
43 56 def compute_cloudpercent(image_path):
44   - array_image, dataset_image = get_raster_as_array(image_path)
45   - cloud = np.sum(array_image == 205)
46   - tot_pix = np.sum(array_image != 254)
47   - return (float(cloud)/float(tot_pix))*100
  57 + array_image, dataset_image = get_raster_as_array(image_path)
  58 + cloud = np.sum(array_image == 205)
  59 + tot_pix = np.sum(array_image != 254)
  60 + return (float(cloud) / float(tot_pix)) * 100
  61 +
48 62  
49 63 def compute_cloud(image):
50   - array, dataset = get_raster_as_array(image)
51   - msk_cloud = (array == 205)
52   - return np.sum(msk_cloud)
  64 + array, dataset = get_raster_as_array(image)
  65 + msk_cloud = (array == 205)
  66 + return np.sum(msk_cloud)
  67 +
53 68  
54 69 def step1(m2_path, m1_path, t0_path, p1_path, p2_path, output_path, ram):
55   - #S(y,x,t) = 1 if (S(y,x,t-1) = 1 and S(y,x,t+1) = 1)
56   - call(["otbcli_BandMath","-ram", str(ram), "-il", m1_path, t0_path, p1_path, "-out", output_path, "-exp", "im2b1==205?((im1b1==100&&im3b1==100)?100:(im1b1==0&&im3b1==0?0:im2b1)):im2b1"])
57   - #S(y,x,t) = 1 if (S(y,x,t-2) = 1 and S(y,x,t+1) = 1)
58   - call(["otbcli_BandMath","-ram", str(ram), "-il", m2_path, output_path, p1_path, "-out", output_path, "-exp", "im2b1==205?((im1b1==100&&im3b1==100)?100:(im1b1==0&&im3b1==0?0:im2b1)):im2b1"])
59   - #S(y,x,t) = 1 if (S(y,x,t-1) = 1 and S(y,x,t+2) = 1)
60   - call(["otbcli_BandMath","-ram", str(ram), "-il", m1_path, output_path, p2_path, "-out", output_path, "-exp", "im2b1==205?((im1b1==100&&im3b1==100)?100:(im1b1==0&&im3b1==0?0:im2b1)):im2b1"])
  70 +
  71 + # Cloud conditions (include pixel flagged as cloud and also as no data
  72 + cloud_nodata_condition = "(im2b1 == 205 || im2b1 == 254)"
  73 +
  74 + # S(y,x,t) = 1 if (S(y,x,t-1) = 1 and S(y,x,t+1) = 1)
  75 + call(
  76 + [
  77 + "otbcli_BandMath",
  78 + "-ram",
  79 + str(ram),
  80 + "-il",
  81 + m1_path,
  82 + t0_path,
  83 + p1_path,
  84 + "-out",
  85 + output_path,
  86 + "-exp",
  87 + cloud_nodata_condition +
  88 + "?((im1b1==100&&im3b1==100)?100:(im1b1==0&&im3b1==0?0:im2b1)):im2b1"])
  89 + # S(y,x,t) = 1 if (S(y,x,t-2) = 1 and S(y,x,t+1) = 1)
  90 + call(
  91 + [
  92 + "otbcli_BandMath",
  93 + "-ram",
  94 + str(ram),
  95 + "-il",
  96 + m2_path,
  97 + output_path,
  98 + p1_path,
  99 + "-out",
  100 + output_path,
  101 + "-exp",
  102 + cloud_nodata_condition +
  103 + "?((im1b1==100&&im3b1==100)?100:(im1b1==0&&im3b1==0?0:im2b1)):im2b1"])
  104 + # S(y,x,t) = 1 if (S(y,x,t-1) = 1 and S(y,x,t+2) = 1)
  105 + call(
  106 + [
  107 + "otbcli_BandMath",
  108 + "-ram",
  109 + str(ram),
  110 + "-il",
  111 + m1_path,
  112 + output_path,
  113 + p2_path,
  114 + "-out",
  115 + output_path,
  116 + "-exp",
  117 + cloud_nodata_condition +
  118 + "?((im1b1==100&&im3b1==100)?100:(im1b1==0&&im3b1==0?0:im2b1)):im2b1"])
  119 +
61 120  
62 121 def step2(t0_path, dem_path, output_path, ram):
63   -
64   - percentage_cloud = compute_cloudpercent(t0_path)
65   - print "cloud percent : " + str(percentage_cloud)
66   - if percentage_cloud < 30:
67   - #S(y,x,t) = 1 if (H(x,y) < Hsmin(t))
68   - hs_min=compute_HSmin(t0_path, dem_path)
69   - print "hs_min: " + str(hs_min)
70   - call(["otbcli_BandMath","-ram", str(ram), "-il", t0_path, dem_path, "-out", output_path, "-exp", "im1b1==205?(im2b1<"+str(hs_min)+"?0:im1b1):im1b1"])
71   - hs_max=compute_HSmax(t0_path, dem_path)
72   - print "hs_max: " + str(hs_max)
73   - #S(y,x,t) = 1 if (H(x,y) > Hsmax(t))
74   - call(["otbcli_BandMath","-ram", str(ram), "-il", output_path, dem_path, "-out", output_path, "-exp", "im1b1==205?(im2b1>"+str(hs_max)+"?100:im1b1):im1b1"])
75   -
  122 +
  123 + percentage_cloud = compute_cloudpercent(t0_path)
  124 + print "cloud percent : " + str(percentage_cloud)
  125 +
  126 + # Perform step 2 only if cloud coverage is less than a threshold value
  127 + # (hard coded for now to 30%)
  128 + cloudpercent_condition = percentage_cloud < 30
  129 +
  130 + if cloudpercent_condition:
  131 + # S(y,x,t) = 1 if (H(x,y) < Hsmin(t))
  132 + hs_min = compute_HSmin(t0_path, dem_path)
  133 + print "hs_min: " + str(hs_min)
  134 + call(["otbcli_BandMath",
  135 + "-ram",
  136 + str(ram),
  137 + "-il",
  138 + t0_path,
  139 + dem_path,
  140 + "-out",
  141 + output_path,
  142 + "-exp",
  143 + "im1b1==205?(im2b1<" + str(hs_min) + "?0:im1b1):im1b1"])
  144 + hs_max = compute_HSmax(t0_path, dem_path)
  145 + print "hs_max: " + str(hs_max)
  146 + # S(y,x,t) = 1 if (H(x,y) > Hsmax(t))
  147 + call(["otbcli_BandMath",
  148 + "-ram",
  149 + str(ram),
  150 + "-il",
  151 + output_path,
  152 + dem_path,
  153 + "-out",
  154 + output_path,
  155 + "-exp",
  156 + "im1b1==205?(im2b1>" + str(hs_max) + "?100:im1b1):im1b1"])
  157 +
  158 + return cloudpercent_condition
  159 +
  160 +
76 161 def step3(t0_path, output_path):
77   -
78   - # four-pixels neighboring
79   - print "Starting step 3"
80   - array, dataset = get_raster_as_array(t0_path)
81   -
82   - #compute 4 pixel snow neighboring
83   - step3_internal(array)
84   -
85   - set_array_as_raster(array, dataset, output_path)
86   - #create file
87   - print "End of step 3"
  162 +
  163 + # four-pixels neighboring
  164 + print "Starting step 3"
  165 + array, dataset = get_raster_as_array(t0_path)
  166 +
  167 + # compute 4 pixel snow neighboring
  168 + step3_internal(array)
  169 +
  170 + set_array_as_raster(array, dataset, output_path)
  171 + # create file
  172 + print "End of step 3"
  173 +
88 174  
89 175 def step3_internal(array):
90   - # Get west, north, east & south elements for [1:-1,1:-1] region of input array
91   - W = array[1:-1,:-2]
92   - N = array[:-2,1:-1]
93   - E = array[1:-1,2:]
94   - S = array[2:,1:-1]
95   -
96   - # Check if all four arrays have 100 for that same element in that region
97   - mask = (W == 100) & (N == 100) & (E == 100) & (S == 100) & (array[1:-1,1:-1] == 205)
98   -
99   - # Use the mask to set corresponding elements
100   - array[1:-1,1:-1][mask] = 100
101   -
  176 + # Get west, north, east & south elements for [1:-1,1:-1] region of input
  177 + # array
  178 + W = array[1:-1, :-2]
  179 + N = array[:-2, 1:-1]
  180 + E = array[1:-1, 2:]
  181 + S = array[2:, 1:-1]
  182 +
  183 + # Check if all four arrays have 100 for that same element in that region
  184 + mask = (W == 100) & (N == 100) & (E == 100) & (
  185 + S == 100) & (array[1:-1, 1:-1] == 205)
  186 +
  187 + # Use the mask to set corresponding elements
  188 + array[1:-1, 1:-1][mask] = 100
  189 +
  190 +
102 191 def step4(t0_path, dem_path, output_path):
103   - # S(y,x,t) = 1 if (S(y+k,x+k,t)(kc(-1,1)) = 1 and H(y+k,x+k)(kc(-1,1)) < H(y,x))
104   - print "Starting step 4"
105   - array, dataset = get_raster_as_array(t0_path)
106   - array_dem, dataset_dem = get_raster_as_array(dem_path)
107   -
108   - # step5
109   - step4_internal(array, array_dem)
110   -
111   - #create file
112   - set_array_as_raster(array, dataset, output_path)
113   - print "End of step 4"
  192 + # S(y,x,t) = 1 if (S(y+k,x+k,t)(kc(-1,1)) = 1 and H(y+k,x+k)(kc(-1,1)) <
  193 + # H(y,x))
  194 + print "Starting step 4"
  195 + array, dataset = get_raster_as_array(t0_path)
  196 + array_dem, dataset_dem = get_raster_as_array(dem_path)
  197 +
  198 + # step5
  199 + step4_internal(array, array_dem)
  200 +
  201 + # create file
  202 + set_array_as_raster(array, dataset, output_path)
  203 + print "End of step 4"
  204 +
114 205  
115 206 def step4_internal(array, array_dem):
116   - # Get 8 neighboring pixels for raster and dem
117   - W = array[1:-1,:-2]
118   - NW = array[:-2,:-2]
119   - N = array[:-2,1:-1]
120   - NE = array[:-2,2:]
121   - E = array[1:-1,2:]
122   - SE = array[2:,2:]
123   - S = array[2:,1:-1]
124   - SW = array[2:,:-2]
125   -
126   - Wdem = array_dem[1:-1,:-2]
127   - NWdem = array_dem[:-2,:-2]
128   - Ndem = array_dem[:-2,1:-1]
129   - NEdem = array_dem[:-2,2:]
130   - Edem = array_dem[1:-1,2:]
131   - SEdem = array_dem[2:,2:]
132   - Sdem = array_dem[2:,1:-1]
133   - SWdem = array_dem[2:,:-2]
134   -
135   - arrdem = array_dem[1:-1,1:-1]
136   - mask =((((W == 100) & (arrdem > Wdem)) | ((N == 100) & (arrdem > Ndem)) | ((E == 100) & (arrdem > Edem)) | ((S == 100) & (arrdem > Sdem)) | ((NW == 100) & (arrdem > NWdem)) | ((NE == 100) & (arrdem > NEdem)) | ((SE == 100) & (arrdem > SEdem)) | ((SW == 100) & (arrdem > SWdem))) & (array[1:-1,1:-1] == 205))
137   -
138   - # Use the mask to set corresponding elements
139   - array[1:-1,1:-1][mask] = 100
  207 + # Get 8 neighboring pixels for raster and dem
  208 + W = array[1:-1, :-2]
  209 + NW = array[:-2, :-2]
  210 + N = array[:-2, 1:-1]
  211 + NE = array[:-2, 2:]
  212 + E = array[1:-1, 2:]
  213 + SE = array[2:, 2:]
  214 + S = array[2:, 1:-1]
  215 + SW = array[2:, :-2]
  216 +
  217 + Wdem = array_dem[1:-1, :-2]
  218 + NWdem = array_dem[:-2, :-2]
  219 + Ndem = array_dem[:-2, 1:-1]
  220 + NEdem = array_dem[:-2, 2:]
  221 + Edem = array_dem[1:-1, 2:]
  222 + SEdem = array_dem[2:, 2:]
  223 + Sdem = array_dem[2:, 1:-1]
  224 + SWdem = array_dem[2:, :-2]
  225 +
  226 + arrdem = array_dem[1:-1, 1:-1]
  227 + mask = (
  228 + (((W == 100) & (
  229 + arrdem > Wdem)) | (
  230 + (N == 100) & (
  231 + arrdem > Ndem)) | (
  232 + (E == 100) & (
  233 + arrdem > Edem)) | (
  234 + (S == 100) & (
  235 + arrdem > Sdem)) | (
  236 + (NW == 100) & (
  237 + arrdem > NWdem)) | (
  238 + (NE == 100) & (
  239 + arrdem > NEdem)) | (
  240 + (SE == 100) & (
  241 + arrdem > SEdem)) | (
  242 + (SW == 100) & (
  243 + arrdem > SWdem))) & (
  244 + array[
  245 + 1:-1, 1:-1] == 205))
  246 +
  247 + # Use the mask to set corresponding elements
  248 + array[1:-1, 1:-1][mask] = 100
  249 +
140 250  
141 251 def compute_stats(image, image_relative, image_reference):
142   - array, dataset = get_raster_as_array(image)
143   - array_relative, dataset = get_raster_as_array(image_relative)
144   - array_reference, dataset = get_raster_as_array(image_reference)
145   - return compute_stats_internal(array, array_relative, array_reference)
146   -
  252 + array, dataset = get_raster_as_array(image)
  253 + array_relative, dataset = get_raster_as_array(image_relative)
  254 + array_reference, dataset = get_raster_as_array(image_reference)
  255 + return compute_stats_internal(array, array_relative, array_reference)
  256 +
  257 +
147 258 def compute_stats_internal(array, array_relative, array_reference):
148   - # Relative Cloud Elimination
149   - msk_cloud_elim = (array_relative == 205) & (array != 205) & (array_reference != 205)
150   - cloud_elim = np.sum(msk_cloud_elim)
151   - # Various stats from paper
152