[GRASS-SVN] r73054 - in grass-addons/grass7/imagery: . i.sentinel.preproc
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Aug 7 06:42:25 PDT 2018
Author: Robifag
Date: 2018-08-07 06:42:24 -0700 (Tue, 07 Aug 2018)
New Revision: 73054
Added:
grass-addons/grass7/imagery/i.sentinel.preproc/
grass-addons/grass7/imagery/i.sentinel.preproc/Makefile
grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.html
grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.py
grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_ES.png
grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_GWF.png
Log:
Added the first version of i.sentinel.preproc module
Added: grass-addons/grass7/imagery/i.sentinel.preproc/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.preproc/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.sentinel.preproc/Makefile 2018-08-07 13:42:24 UTC (rev 73054)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.sentinel.preproc
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.html
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.html (rev 0)
+++ grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.html 2018-08-07 13:42:24 UTC (rev 73054)
@@ -0,0 +1,324 @@
+<h2>DESCRIPTION</h2> <em>i.sentinel.preproc</em> allows to import Sentinel 2
+images and perform atmospheric correction. <p> i.sentinel.preproc is a module
+for the preprocessing of Sentinel 2 images (Level-1C Single Tile product) which
+wraps the import and the atmospheric correction using respectively <a
+href="i.sentinel.import.html">i.sentinel.import</a> and <a
+href="i.atcorr.html">i.atcorr</a>.<br>
+It works both with Sentinel 2A and 2B images.<br>
+The aim is to provide a simplified module which allows importing images, which
+area downloaded using <a
+href="i.sentinel.download.html">i.sentinel.download</a> or any other sources,
+and performing the atmospheric correction avoiding users to provide all the
+required input parameters manually. In fact, regarding the atmospheric
+correction performed with i.atcorr one of the most challenging steps,
+especially for unexperienced users, is the compiling of the control file with
+all the required parameters to parametrize the 6S (<em>Second Simulation of
+Satellite Signal in the Solar Spectrum</em>) model on which i.atcorr is
+based.<br><br> To run i.atcorr, users have to provide the so-called control
+file in which all the parameters (geometrical conditions, date, time, longitude
+and latitude of the center of the scene, atmospheric model, aerosol model,
+visibility or Aerosol Optical Depth -AOD- value, mean elevation target and
+bands number) have to be specified with precise syntax rules and codes.<br>
+i.sentinel.preproc retrieves as many parameters as possible from the metadata
+file (e.g. Geometrical conditions, data and time and bands number), longitude
+and latitude are automatically computed from the computational region while
+others like the mean target elevation above sea level from the input digital
+elevation model (DEM). Only a few parameters have to be provided by users who
+can choose the proper option from a drop-down menu thus avoiding to enter the
+corresponding code. In any case, i.sentinel.preproc writes a temporary control
+file, changing it according to the band number, following the syntax rules and
+codes of i.atcorr and then it runs i.atcorr for all bands.<br><br>
+<br><br>
+<center>
+<img src="i_sentinel_preproc_GWF.png" width="30%">
+<br><br>
+<i>Fig: Module General WorkFlow</i>
+</center>
+<br><br>
+When all bands have been processed by the integrated version of i.atcorr, an
+histogram equalization grayscale color scheme is applied.
+<br><br>
+If the <b>-t</b> flag is checked, a text file ready to be used as input for <a
+href="i.sentinel.mask">i.sentinel.mask</a> will be created. In this case a name
+for the output text file has to be specified.
+<br><br>
+NOTE: as for i.atcorr, current region settings are ignored. The region is
+temporary set to maximum image extent and restored at the end of the process.
+<br><br>
+<em><b>Important</b></em>: i.sentinel.preproc requires all the bands of a
+Sentinel 2 images. If the module is used only for the atmospheric correction,
+all bands from *_B01 to *_B12 must be imported.<br>Moreover, the original bands
+name has to be kept unchanged (e.g if the original name is
+<em>T17SPV_20180315T160021_B02</em> the imported raster map in the GIS DATABSE
+must be named <em>T17SPV_20180315T160021_B02</em>).
+<br><br>
+<h3>Import</h3>
+<p>
+i.sentinel.preproc allows the import of all the bands of a Sentinel 2 image.
+The required input is the <b>.SAFE folder</b> downloaded using
+i.sentinel.download or any other source (e.g. Copernicus Open Access Hub). Note
+that in the case that spatial reference system of input data differs from GRASS
+location, the input data are reprojected.<br> Number of imported bands <b>can
+not</b> be reduced, all bands are automatically imported by default.<br><br>
+
+<em><b>Important</b></em>: i.sentinel.preproc allows the import of one image at
+a time because the input .SAFE folder is also used to automatically identify
+the corresponding metadata file that is used during the atmospheric
+correction.<br><br>
+
+The import can be skipped using the <b>-i</b> flag. Note that even if the
+import is skipped the input .SAFE folder must be specified to automatically
+retrieve the metadata file.
+
+<h3>Atmospheric correction</h3>
+<p>
+i.sentinel.preproc allows performing atmospheric correction of all bands of a
+Sentinel 2 scene with a single process using i.atcorr. Unlike i.atcorr, it
+writes the control file changing it according to the band number. The only
+required inputs are:
+<ul>
+ <li><b>input_dir</b> = the *.SAFE directory where the image and metadata file (MTD_MSIL1C.xml or S2A_OPER_MTD_SAFL1C_PDMC_*.xml depending on naming convention) are stored,
+ <li><b>elevation</b> = raster of a digital elevation model,
+ <li><b>visibility or AOD value</b> = raster of a visibility map or an AOD value (<em>see AOD section</em>),
+ <li><b>Atmospheric model</b> = to be choosen from the drop-down menu,
+ <li><b>Aerosol model</b> = to be choosen from the drop-down menu,
+ <li><b>suffix</b> = a suffix for the output maps name
+ <li><b>rescale</b> = the output range of values for the corrected bands, for example: 0-255, 0-1, 1-10000 (default value 0-1).
+</ul>
+<br>
+The module writes the control file automatically starting from the input above.
+
+<h4>Control file</h4>
+<p>
+i.atcorr requires a control file to parametrize the 6S algorithm on which it is
+based.
+
+<p>
+Below an example of the control file, taken from the i.atcorr manual page, of a
+Sentinel 2A image:
+
+<div class="code"><pre>
+25 - geometrical conditions = Sentinel 2A
+5 4 19.737 -78.727 35.748 - month day hh.ddd longitude latitude ("hh.ddd" is in decimal hours GMT)
+2 - atmospheric model = midlatitude summer
+1 - aerosols model = continental
+0 - visibility [km] (aerosol model concentration)
+0.07 - AOD at 550nm
+-0.124 - mean target elevation above sea level [km]
+-1000 - sensor height (here, sensor on board a satellite)
+167 - sensor band = Sentinel2A Blue band B2
+</pre></div>
+
+Using i.sentinel.preproc the only parameters from the list above that users
+have to provide are: atmospheric model, aerosol model, visibility or AOD value.
+The others are automatically retrieved from the metadata file, input elevation
+map and bands.
+<br><br>
+<ol>
+<li><b>Geometrical conditions</b>
+<p>
+The geometrical condition of the satellite are read from the metadata file and
+converted to the corresponding i.atcorr code, 25 for Sentinel 2A mission and 26
+for Sentinel 2B.
+<br><br>
+<li><b> Date, time, longitude and latitude </b>
+<p>
+Date (month and day) and time are read from the metadata file. The date (with
+the format YYYY-MM-DDTHH:MM:SSZ) is converted in a standard format and only the
+month and the day are selected and added to the control file.<br><br> Time is
+already in Greenwich Mean Time (GMT), as i.atcorr requires, and it's
+automatically converted to decimal hours.<br> Longitude and latitude are
+computed from the computational region and converted to WGS84 decimal
+coordinates.
+<br><br>
+<li><b> Atmospheric model</b>
+<p>
+Only some options are available:
+<ul>
+ <li>Automatic
+ <li>No gaseous absorption
+ <li>Tropical
+ <li>Midlatitude summer
+ <li>Midlatitude winter
+ <li>Subarctic summer
+ <li>Subarctic winter
+ <li>Us standard 62
+</ul><br>
+Users can choose the proper option from a drop-down menu. The desired model is
+automatically converted to the corresponding code and added to the control
+file.
+<br><br><em><b>Automatic</b> option</em><br>
+The default option is <em>Automatic</em> which consists in the automatic
+identification of the proper atmospheric model for the input image. The
+<em>Automatic</em> option reads the latitude of the center of the computational
+region and uses it to choose between Midlatitude (15.00 > lat <= 45.00),
+Tropical (-15.00 > lat <= 15.00) and Subarctic (45.00 > lat <= 60.00) for
+Northern Hemisphere (obviously it also works for the Southern Hemisphere).
+Then, the month from the acquisition date is used to distinguish summer or
+winter in case of Midlatitude or Subarctic model. Once the proper atmospheric
+model is identified, it is converted to the corresponding code and added to the
+control file.<br> Note that this is a simplified and standardized method to
+identify the atmospheric model. Obviously, it is possible to choose other
+options from those available.
+<br><br>
+<li><b> Aerosol model</b><br>
+<p>
+Also in this case, only some options are available and users have to select the
+desired one from the drop-down menu, then it is converted to the corresponding
+code and added to the control file.
+<ul>
+ <li>no aerosols
+ <li>continental model
+ <li>maritime model
+ <li>urban model
+ <li>shettle model for background desert aerosol
+ <li>biomass burning
+ <li>stratospheric model
+</ul><br>
+No automatic procedure has been implemented in this case.
+<br><br>
+<li><b> Visibility or AOD </b><br>
+<p>
+By default, i.sentinel.preproc uses the input visibility map to estimate a
+visibility value to be added in the control file. If no visibility map is
+available for the processed scene, it is possible to use an estimated Aerosol
+Optical Depth (AOD) value checking the <b>-a</b> flag.<br>
+If the -a flag is checked and a visibility map is provided, the visibility will
+be ignored and no mean visibility value will be computed and added to the
+control file. Whereas, if the -a flag isn't checked and an AOD value is
+provided it will be ignored and not added to the control file.<br><br> In the
+same way, if the -a flag is checked and a visibility map is provided it will be
+excluded from atmospheric correction process.<br><br>
+<b> AOD</b><br><br>
+The AOD value can be specified by users (e.g. <tt>aod_value=0.07</tt>) or
+automatically retrieved from an AERONET file to be given as input instead of
+the AOD value.<br> i.sentinel.preproc reads the AERONET file, identify the
+closest available date to the scene date and compute AOD at 550nm using the
+closest upper and lower wavelength to 550 (e.g. 500nm and 675nm) and applying
+the Angstrom coefficient. <br><br>
+The type of AERONET file is a Combined file for All Points (Level 1.5 or
+2.0)<br> To download this kind of file:<br>
+<ol>
+<li>Go to <a href=http://aeronet.gsfc.nasa.gov/cgi-bin/webtool_opera_v2_inv>http://aeronet.gsfc.nasa.gov</a>
+<li>Choose the site you want to get data from
+<li>Choose the data you want to get data for
+<li>Tick the box near the bottom labelled as 'Combined file (all products without phase functions)'
+<li>Choose either Level 1.5 or Level 2.0 data. Level 1.5 data is unscreened, so contains far more data meaning it is more likely for users to find data near your specified time
+<li>Choose 'All Points' under Data Format
+<li>Download the file
+<li>Unzip (the file has a .dubovik extension)
+</ol>
+<br><br>
+Then, giving this file as input (e.g.
+<tt>aeronet_file=your_path/*.dubovik</tt>), the AOD at 550nm will be
+automatically computed and added to the control file.<br><br> NOTE: as i.atcorr
+manual explain, if an AOD value is provided a value 0 for the visibility has to
+be entered with the AOD value in the following line. Obviously,
+i.sentinel.preproc takes into account this syntax rule and automatically adds a
+0 value for visibility (or -1 if AOD=0) if an AOD value is provided (through
+both <tt>aod_value</tt> and <tt>aeronet_file</tt>).
+<br><br>
+<li><b> Mean target elevation above sea level </b><br>
+<p>
+Mean target elevation above sea level is automatically estimated from the input
+digital elevation model. According to the rules for writing the contol file of
+i.atcorr, the mean elevation value is added as a negative value and converted
+in kilometers (e.g. if mean=121 in the control file it will be written in
+[-km], i.e., -0.121).
+<br><br>
+<li><b> Sensor height </b><br>
+<p>
+Since the sensor is on board a satellite, the sensor height is automatically
+set to -1000.
+<br><br>
+<li><b> Sensor band </b><br>
+<p>
+The number of the band changes automatically according to the band that is
+processed at that moment. The module reads the name of the band and converts it
+into the corresponding code.
+</ol>
+
+<h2>EXAMPLE</h2>
+<p>
+The example illustrates how to run i.sentinel.preproc for a Sentinel 2A image
+(S2A_MSIL1C_20180315T160021_N0206_R097_T17SPV_20180315T194425.SAFE) in the
+North Carolina location.<br> The AERONET file has been downloaded from the
+<em>EPA-Res_Triangle_Pk</em> station.
+<div class="code"><pre>
+i.sentinel.preproc -a -t input_dir=/path/S2A_MSIL1C_20180315T160021_N0206_R097_T17SPV_20180315T194425.SAFE elevation=elevation atmospheric_model=Automatic aerosol_model="Continental model" aeronet_file=path/180301_180331_EPA-Res_Triangle_Pk.dubovik suffix=cor text_file=/path/input_cloud_mask.txt
+<pre></div>
+<p>
+Here is the control file automatically written for Band 02 of the input scene
+<div class="code"><pre>
+25
+5 4 19.74 -78.728 35.749
+3 -The Automatic option identified the Midlatitude Winter as the proper model for the scene
+1
+0 -The visibility is set to 0 with AOD in the following line
+0.18867992317 -AOD computed from the input AERONET file
+-0.122
+-1000
+167
+</pre></div>
+Here is the output text file ready to be used as input for i.sentinel.mask (-t
+flag)
+<div class="code"><pre>
+blue=T17SPV_20180315T160021_B02_cor
+green=T17SPV_20180315T160021_B03_cor
+red=T17SPV_20180315T160021_B04_cor
+swir11=T17SPV_20180315T160021_B11_cor
+nir=T17SPV_20180315T160021_B08_cor
+swir12=T17SPV_20180315T160021_B12_cor
+nir8a=T17SPV_20180315T160021_B8A_cor
+</pre></div>
+<br>
+<div align="center" style="margin: 10px">
+<a href="i_sentinel_preproc_ES.png">
+<img src="i_sentinel_preproc_ES.png" width="538" height="444" alt="i.atcorr example" border="0">
+</a><br><br>
+<i>Figure: Sentinel-2A Band 02</i>
+</div><br>
+
+<h2>REQUIREMENTS</h2>
+<ul>
+<li><a href="i.sentinel.import.html">i.sentinel.import</a>
+</ul>
+
+<h2>IMPORTANT NOTES</h2>
+<ul>
+<li>i.sentinel.preproc integrates a simplyfied version of both modules
+(i.sentinel.import and i.atcorr), only some options are available. For
+instance, if it's necessary a strong customization (e.g. definition of your own
+atmospheric or aerosol model), please refer to i.atcorr.<br><br>
+<li>i.sentinel.preproc works with Sentinel 2 images whose names follow both the
+New Compact Naming Convention (e.g.
+S2A_MSIL1C_20170527T102031_N0205_R065_T32TMQ_20170527T102301.SAFE) and the Old
+Format Naming Convention (e.g.
+S2A_OPER_PRD_MSIL1C_PDMC_20160930T155112_R079_V20160930T095022_20160930T095944.SAFE).
+For further information about the naming convention see <a
+href="https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/naming-convention">ESA
+Sentinel User Guide</a>.
+</ul>
+
+<h2>FOLLOW UP</h2>
+<ul>
+<li>Implement download functionality avoiding dependencies
+<li>Integrate the Topographic Correction
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.sentinel.download.html">i.sentinel.download</a>,
+<a href="i.sentinel.import.html">i.sentinel.import</a>,
+<a href="i.atcorr.html">i.atcorr</a>,
+<a href="i.sentinel.import.html">i.sentinel.mask</a>,
+<a href="r.import.html">r.import</a>,
+<a href="r.extenal.html">r.external</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Roberta Fagandini, GSoC 2018 student<br>
+<a href="https://wiki.osgeo.org/wiki/User:Mlennert">Moritz Lennert</a><br>
+<a href="https://wiki.osgeo.org/wiki/User:Robertomarzocchi">Roberto Marzocchi</a>
Added: grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.py
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.py (rev 0)
+++ grass-addons/grass7/imagery/i.sentinel.preproc/i.sentinel.preproc.py 2018-08-07 13:42:24 UTC (rev 73054)
@@ -0,0 +1,706 @@
+#!/usr/bin/env python
+# coding=utf-8
+#
+############################################################################
+#
+# MODULE: i.sentinel.preproc
+# AUTHOR(S): Roberta Fagandini, Moritz Lennert, Roberto Marzocchi
+# PURPOSE: Import and perform atmospheric correction for Sentinel-2 images
+#
+# COPYRIGHT: (C) 2018 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+############################################################################
+
+#%Module
+#% description: Import and perform atmospheric correction for Sentinel-2 images
+#% overwrite: yes
+#% keywords: imagery
+#% keywords: sentinel
+#% keywords: download
+#% keywords: import
+#% keywords: atmospheric correction
+#%End
+#%option
+#% key: input_dir
+#% type: string
+#% gisprompt: old,dir,dir
+#% description: Name of the directory where the image and metadata file are stored (*.SAFE)
+#% required : yes
+#% multiple: no
+#% guisection: Input
+#%end
+#%option
+#% key: elevation
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of input elevation raster map (in m)
+#% required : yes
+#% multiple: no
+#% guisection: Input
+#%end
+#%option
+#% key: visibility
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of input visibility raster map (in Km)
+#% required : no
+#% multiple: no
+#% guisection: Input
+#%end
+#%option
+#% key: atmospheric_model
+#% type: string
+#% description: Select the proper Atmospheric model
+#% options: Automatic,No gaseous absorption,Tropical,Midlatitude summer,Midlatitude winter,Subarctic summer,Subarctic winter,Us standard 62
+#% answer: Automatic
+#% required : yes
+#% multiple: no
+#% guisection: 6S Parameters
+#%end
+#%option
+#% key: aerosol_model
+#% type: string
+#% description: Select the proper Aerosol model
+#% options: No aerosols,Continental model,Maritime model,Urban model,Shettle model for background desert aerosol,Biomass burning,Stratospheric model
+#% answer: Continental model
+#% required : yes
+#% multiple: no
+#% guisection: 6S Parameters
+#%end
+#%option
+#% key: aod_value
+#% type: string
+#% description: AOD value at 550nm
+#% required : no
+#% guisection: 6S Parameters
+#%end
+#%option
+#% key: aeronet_file
+#% type: string
+#% gisprompt: old,file,file
+#% description: Name of the AERONET file for computing AOD at 550nm
+#% required : no
+#% multiple: no
+#% guisection: 6S Parameters
+#%end
+#%option
+#% key: suffix
+#% type: string
+#% description: Suffix for output maps
+#% required : yes
+#% guisection: Output
+#%end
+#%option
+#% key: rescale
+#% key_desc: min,max
+#% type: string
+#% description: Rescale output raster map
+#% answer: 0,1
+#% required : no
+#% guisection: Output
+#%end
+#%option
+#% key: text_file
+#% type: string
+#% gisprompt: new,file,file
+#% description: Name of of output text file to be used as input in i.sentinel.mask
+#% required : no
+#% guisection: Output
+#%end
+#%flag
+#% key: a
+#% description: Use AOD instead visibility
+#% guisection: 6S Parameters
+#%end
+#%flag
+#% key: t
+#% description: Create the input text file for i.sentinel.mask
+#% guisection: Output
+#%end
+#%flag
+#% key: r
+#% description: Reproject raster data using r.import if needed
+#% guisection: Input
+#%end
+#%flag
+#% key: i
+#% description: Skip import of Sentinel bands
+#% guisection: Input
+#%end
+
+import grass.script as gscript
+import xml.etree.ElementTree as et
+from datetime import datetime
+import os
+import math
+import sys
+import shutil
+import re
+import glob
+import atexit
+
+
+def main ():
+
+
+ bands = {}
+ cor_bands = {}
+ dem = options['elevation']
+ vis = options['visibility']
+ input_dir = options['input_dir']
+ check_ndir = 0
+ check_odir = 0
+ # Check if the input folder has old or new name
+ # Check if the input folder belongs to a L1C image
+ level_dir = os.path.basename(input_dir).split('_')
+ if level_dir[1] == 'OPER' and level_dir[3] == 'MSIL1C':
+ check_odir = 1
+ filename = [i for i in os.listdir(input_dir) if i.startswith("S")]
+ string = str(filename).strip("['']")
+ mtd_file = options['input_dir'] + '/' + string
+ elif level_dir[1] == 'MSIL1C':
+ check_ndir = 1
+ mtd_file = options['input_dir'] + '/MTD_MSIL1C.xml'
+ else:
+ gscript.fatal('The input directory does not belong to a L1C Sentinel image. Please check the input directory')
+ # Check if Metadata file exists
+ if not os.path.isfile(mtd_file):
+ gscript.fatal('Metadata file not found. Please check the input directory')
+ atmo_mod = options['atmospheric_model']
+ aerosol_mod = options['aerosol_model']
+ aeronet_file = options['aeronet_file']
+ check_file = 0
+ check_value = 0
+ mapset = gscript.gisenv()['MAPSET']
+ suffix = options['suffix']
+ rescale = options['rescale']
+ processid = os.getpid()
+ txt_file = options['text_file']
+ tmp_file = gscript.tempfile()
+
+ # Import bands
+ if not flags["i"]:
+ try:
+ if flags["r"]:
+ gscript.run_command('i.sentinel.import',
+ input=options['input_dir'],
+ flags='r')
+ else:
+ gscript.run_command('i.sentinel.import',
+ input=options['input_dir'])
+ except:
+ gscript.fatal('Module rquire i.sentinel.import. Please install it using g.extension.')
+
+ # Create xml "tree" for reading parameters from metadata
+ tree = et.parse(mtd_file)
+ root = tree.getroot()
+
+ # Start reading the xml file
+ if check_ndir == 1:
+ for elem in root[0].findall('Product_Info'):
+ datatake = elem.find('Datatake')
+ # Geometrical conditions = sensor
+ sensor = datatake.find('SPACECRAFT_NAME')
+ # Acquisition date and time
+ time_str = elem.find('GENERATION_TIME')
+ # Date and time conversion
+ time_py = datetime.strptime(time_str.text,'%Y-%m-%dT%H:%M:%S.%fZ')
+ # Compute decimal hour
+ dec_hour = (float(time_py.hour)
+ + float(time_py.minute)/60
+ + float(time_py.second)/3600)
+ # Read input bands from metadata
+ product = elem.find('Product_Organisation')
+ g_list = product.find('Granule_List')
+ granule = g_list.find('Granule')
+ images = granule.find('IMAGE_FILE')
+ img_name = images.text.split('/')
+ # Check if input exist and if the mtd file corresponds with the input image
+ for img in root.iter('IMAGE_FILE'):
+ a = img.text.split('/')
+ b = a[3].split('_')
+ if gscript.find_file(a[3],
+ element = 'cell',
+ mapset = mapset)['file'] or b[2] == 'TCI':
+ if b[2] == 'B01':
+ bands['costal'] = a[3]
+ elif b[2] == 'B02':
+ bands['blue'] = a[3]
+ elif b[2] == 'B03':
+ bands['green'] = a[3]
+ elif b[2] == 'B04':
+ bands['red'] = a[3]
+ elif b[2] == 'B05':
+ bands['re5'] = a[3]
+ elif b[2] == 'B06':
+ bands['re6'] = a[3]
+ elif b[2] == 'B07':
+ bands['re7'] = a[3]
+ elif b[2] == 'B08':
+ bands['nir'] = a[3]
+ elif b[2] == 'B8A':
+ bands['nir8a'] = a[3]
+ elif b[2] == 'B09':
+ bands['vapour'] = a[3]
+ elif b[2] == 'B10':
+ bands['cirrus'] = a[3]
+ elif b[2] == 'B11':
+ bands['swir11'] = a[3]
+ elif b[2] == 'B12':
+ bands['swir12'] = a[3]
+ else:
+ gscript.fatal(("One or more input bands are missing or \n the metadata file belongs to another image ({}).").format(
+ img_name[3].replace('_B01','')))
+
+ if check_odir == 1:
+ for elem in root[0].findall('Product_Info'):
+ datatake = elem.find('Datatake')
+ # Geometrical conditions = sensor
+ sensor = datatake.find('SPACECRAFT_NAME')
+ # Acquisition date and time
+ time_str = elem.find('GENERATION_TIME')
+ # Date and time conversion
+ time_py = datetime.strptime(time_str.text,'%Y-%m-%dT%H:%M:%S.%fZ')
+ # Compute decimal hour
+ dec_hour = (float(time_py.hour)
+ + float(time_py.minute)/60
+ + float(time_py.second)/3600)
+ # Read input bands from metadata
+ product = elem.find('Product_Organisation')
+ g_list = product.find('Granule_List')
+ granule = g_list.find('Granules')
+ images = granule.find('IMAGE_ID')
+ # Check if input exist and if the mtd file corresponds with the input image
+ for img in root.iter('IMAGE_ID'):
+ b = img.text.split('_')
+ if gscript.find_file(img.text,
+ element = 'cell',
+ mapset = mapset)['file']:
+ if b[10] == 'B01':
+ bands['costal'] = img.text
+ elif b[10] == 'B02':
+ bands['blue'] = img.text
+ elif b[10] == 'B03':
+ bands['green'] = img.text
+ elif b[10] == 'B04':
+ bands['red'] = img.text
+ elif b[10] == 'B05':
+ bands['re5'] = img.text
+ elif b[10] == 'B06':
+ bands['re6'] = img.text
+ elif b[10] == 'B07':
+ bands['re7'] = img.text
+ elif b[10] == 'B08':
+ bands['nir'] = img.text
+ elif b[10] == 'B8A':
+ bands['nir8a'] = img.text
+ elif b[10] == 'B09':
+ bands['vapour'] = img.text
+ elif b[10] == 'B10':
+ bands['cirrus'] = img.text
+ elif b[10] == 'B11':
+ bands['swir11'] = img.text
+ elif b[10] == 'B12':
+ bands['swir12'] = img.text
+ else:
+ gscript.fatal(("One or more input bands are missing or \n the metadata file belongs to another image ({}).").format(
+ images.text.replace('_B09','')))
+
+ # Check if input exist
+ for key, value in bands.items():
+ if not gscript.find_file(value,
+ element = 'cell',
+ mapset = mapset)['file']:
+ gscript.fatal(('Raster map <{}> not found.').format(value))
+
+ # Check if output already exist
+ for key, value in bands.items():
+ if not os.getenv('GRASS_OVERWRITE'):
+ if gscript.find_file(value + '_' + suffix,
+ element = 'cell',
+ mapset = mapset)['file']:
+ gscript.fatal(('Raster map {} already exists.').format(
+ value + '_' + suffix))
+
+ # Check if output name for the text file has been specified
+ if flags["t"]:
+ if options['text_file'] == '':
+ gscript.fatal('Output name is required for the text file. Please specified it')
+
+ # Set temp region to image max extent
+ gscript.use_temp_region()
+ gscript.run_command('g.region',
+ rast=bands.values(),
+ flags='a')
+ gscript.message(_("--- The computational region has been temporarily set to image max extent ---"))
+
+ if flags["a"]:
+ if vis!='':
+ if options['aod_value']!='' and aeronet_file!='':
+ gscript.warning(_('--- Visibility map will be ignored ---'))
+ gscript.fatal('Only one parameter must be provided, AOD value or AERONET file')
+ elif options['aod_value']=='' and aeronet_file=='':
+ gscript.warning(_('--- Visibility map will be ignored ---'))
+ gscript.fatal('If -a flag is checked an AOD value or AERONET file must be provided')
+ elif options['aod_value']!='':
+ gscript.warning(_('--- Visibility map will be ignored ---'))
+ check_value = 1
+ aot550 = options['aod_value']
+ elif aeronet_file!='':
+ gscript.warning(_('--- Visibility map will be ignored ---'))
+ elif options['aod_value']!='' and aeronet_file!='':
+ gscript.fatal('Only one parameter must be provided, AOD value or AERONET file')
+ elif options['aod_value']!='':
+ check_value = 1
+ aot550 = options['aod_value']
+ elif aeronet_file!='':
+ gscript.message(_('--- Computing AOD from input AERONET file ---'))
+ elif options['aod_value']=='' and aeronet_file=='':
+ gscript.fatal('If -a flag is checked an AOD value or AERONET file must be provided')
+ else:
+ if vis!='':
+ if options['aod_value']!='' or aeronet_file!='':
+ gscript.warning(_('--- AOD will be ignored ---'))
+ check_file = 1
+ stats_v = gscript.parse_command('r.univar', flags='g', map=vis)
+ vis_mean = int(float(stats_v['mean']))
+ gscript.message(_('--- Computed visibility mean value: {} Km ---'.format(vis_mean)))
+ elif vis=='' and (options['aod_value']!='' or aeronet_file!=''):
+ gscript.fatal('Check the -a flag to use AOD instead of visibility')
+ else:
+ gscript.fatal('No visibility map has been provided')
+
+ # Retrieve longitude and latitude of the centre of the computational region
+ c_region = gscript.parse_command('g.region', flags='bg')
+ lon = (float(c_region['ll_clon']))
+ lat = (float(c_region['ll_clat']))
+
+ # Read and compute AOD from AERONET file
+ if check_value == 0 and check_file == 0:
+ i=0
+ cc=0
+ count=0
+ columns = []
+ m_time = []
+ dates_list = []
+ t_columns = []
+ i_col = []
+ coll = []
+ wl = []
+ for row in file(aeronet_file):
+ count+=1
+ if count==4:
+ columns = row.split(',')
+ # Search for the closest date and time to the acquisition one
+ count=0
+ for row in file(aeronet_file):
+ count+=1
+ if count>=5:
+ columns = row.split(',')
+ m_time.append(columns[0] + ' ' + columns[1])
+
+ dates = [datetime.strptime(row, '%d:%m:%Y %H:%M:%S') for row in m_time]
+ dates_list.append(dates)
+ format_bd = time_py.strftime('%d/%m/%Y %H:%M:%S')
+ base_date = str(format_bd)
+ b_d = datetime.strptime(base_date,'%d/%m/%Y %H:%M:%S')
+
+ for line in dates_list:
+ closest = min(line, key=lambda x: abs(x - b_d))
+ timedelta = abs(closest - b_d)
+ # Search for the closest wavelengths (upper and lower) to 550
+ count=0
+ for row in file(aeronet_file):
+ count+=1
+ if count==4:
+ t_columns = row.split(',')
+ for i, col in enumerate(t_columns):
+ if "AOT_" in col:
+ i_col.append(i)
+ coll.append(col)
+ for line in coll:
+ l = line.split('_')
+ wl.append(int(l[1]))
+
+ aot_req = 550
+ upper = min([ i for i in wl if i >= aot_req], key=lambda x:abs(x-aot_req))
+ lower = min([ i for i in wl if i < aot_req], key=lambda x:abs(x-aot_req))
+
+ count=0
+ for row in file(aeronet_file):
+ count+=1
+ if count==dates.index(closest)+5:
+ t_columns = row.split(',')
+ count2=0
+ check_up=0
+ check_lo=0
+ while count2<len(i_col) and check_up<1:
+ # Search for the not null value for the upper wavelength
+ if t_columns[wl.index(upper)+i_col[0]]=="N/A":
+ aot_req_tmp = upper
+ upper = min([ i for i in wl if i > aot_req_tmp], key=lambda x:abs(x-aot_req_tmp))
+ else:
+ wl_upper = float(upper)
+ aot_upper = float(t_columns[wl.index(upper)+i_col[0]])
+ check_up=1
+ count2+=1
+ count2=0
+ while count2<len(i_col) and check_lo<1:
+ # Search for the not null value for the lower wavelength
+ if t_columns[wl.index(lower)+i_col[0]]=="N/A":
+ aot_req_tmp = lower
+ lower = min([ i for i in wl if i < aot_req_tmp], key=lambda x:abs(x-aot_req_tmp))
+ else:
+ wl_lower = float(lower)
+ aot_lower = float(t_columns[wl.index(lower)+i_col[0]])
+ check_lo=1
+ count2+=1
+ # Compute AOD at 550 nm
+ alpha = math.log(aot_lower/aot_upper)/math.log(wl_upper/wl_lower)
+ aot550 = math.exp(math.log(aot_lower) - math.log(550.0/wl_lower)*alpha)
+ gscript.message(_('--- Computed AOD at 550 nm: {} ---'.format(aot550)))
+
+ # Compute mean target elevation in km
+ stats_d = gscript.parse_command('r.univar', flags='g', map=dem)
+ mean = (float(stats_d['mean']))
+ conv_fac = -0.001
+ dem_mean = mean * conv_fac
+ gscript.message(_('--- Computed mean target elevation above sea level: {:.3f} m ---'.format(mean)))
+
+ # Start compiling the control file
+ for key, bb in bands.items():
+ gscript.message(_('--- Compiling the control file.. ---'))
+ text = open(tmp_file, "w")
+ # Geometrical conditions
+ if sensor.text == 'Sentinel-2A':
+ text.write(str(25) + "\n")
+ elif sensor.text == 'Sentinel-2B':
+ text.write(str(26) + "\n")
+ else:
+ gscript.fatal('The input image does not seem to be a Sentinel image')
+ text.write('{} {} {:.2f} {:.3f} {:.3f}'.format(
+ time_py.month,
+ time_py.day,
+ dec_hour,
+ lon,
+ lat) + "\n")
+ # Atmospheric model
+ winter = [1, 2, 3, 4, 10, 11, 12]
+ summer = [5, 6, 7, 8, 9]
+ if atmo_mod == 'Automatic':
+ if lat > -15.00 and lat <= 15.00: # Tropical
+ text.write('1' + "\n")
+ elif lat > 15.00 and lat <= 45.00:
+ if time_py.month in winter: # Midlatitude winter
+ text.write('3' + "\n")
+ else: # Midlatitude summer
+ text.write('2' + "\n")
+ elif lat > -15.00 and lat <= -45.00:
+ if time_py.month in winter: # Midlatitude summer
+ text.write('2' + "\n")
+ else: # Midlatitude winter
+ text.write('3' + "\n")
+ elif lat > 45.00 and lat <= 60.00:
+ if time_py.month in winter: # Subarctic winter
+ text.write('5' + "\n")
+ else: # Subartic summer
+ text.write('4' + "\n")
+ elif lat > -45.00 and lat <= -60.00:
+ if time_py.month in winter: # Subarctic summer
+ text.write('4' + "\n")
+ else: # Subartic winter
+ text.write('5' + "\n")
+ else:
+ gscript.fatal('Latitude {} is out of range'.format(lat))
+ elif atmo_mod == 'No gaseous absorption':
+ text.write('0' + "\n") # No gas abs model
+ elif atmo_mod == 'Tropical':
+ text.write('1' + "\n") # Tropical model
+ elif atmo_mod == 'Midlatitude summer':
+ text.write('2' + "\n") # Mid sum model
+ elif atmo_mod == 'Midlatitude winter':
+ text.write('3' + "\n") # Mid win model
+ elif atmo_mod == 'Subarctic summer':
+ text.write('4' + "\n") # Sub sum model
+ elif atmo_mod == 'Subarctic winter':
+ text.write('5' + "\n") # Sub win model
+ elif atmo_mod == 'Us standard 62':
+ text.write('6' + "\n") # Us 62 model
+ # Aerosol model
+ if aerosol_mod == 'No aerosols':
+ text.write('0' + "\n") # No aerosol model
+ elif aerosol_mod == 'Continental model':
+ text.write('1' + "\n") # Continental aerosol model
+ elif aerosol_mod == 'Maritime model':
+ text.write('2' + "\n") # Maritimw aerosol model
+ elif aerosol_mod == 'Urban model':
+ text.write('3' + "\n") # Urban aerosol model
+ elif aerosol_mod == 'Shettle model for background desert aerosol':
+ text.write('4' + "\n") # Shettle aerosol model
+ elif aerosol_mod == 'Biomass burning':
+ text.write('5' + "\n") #Biomass aerosol model
+ elif aerosol_mod == 'Stratospheric model':
+ text.write('6' + "\n") # Stratospheric aerosol model
+ # Visibility and/or AOD
+ if not flags["a"] and vis != '':
+ text.write('{}'.format(vis_mean) + "\n")
+ elif flags["a"] and vis != '':
+ if aot550 != 0:
+ text.write('0' + "\n") # Visibility
+ text.write('{}'.format(aot550) + "\n")
+ elif aot550 == 0:
+ text.write('-1' + "\n") # Visibility
+ text.write('{}'.format(aot550) + "\n")
+ elif vis == '' and aot550 != 0:
+ text.write('0' + "\n") # Visibility
+ text.write('{}'.format(aot550) + "\n")
+ elif vis == '' and aot550 == 0:
+ text.write('-1' + "\n") # Visibility
+ text.write('{}'.format(aot550) + "\n")
+ else:
+ gscript.fatal('Unable to retrieve visibility or AOD value, check the input')
+ text.write('{:.3f}'.format(dem_mean) + "\n") # Mean elevation
+ text.write('-1000' + "\n") # Sensor height
+ # Band number
+ b = bb.split('_')
+ if check_ndir == 1:
+ band_n = b[2]
+ else:
+ band_n = b[10]
+ if band_n == 'B01' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('166')
+ elif band_n == 'B02' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('167')
+ elif band_n == 'B03' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('168')
+ elif band_n == 'B04' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('169')
+ elif band_n == 'B05' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('170')
+ elif band_n == 'B06' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('171')
+ elif band_n == 'B07' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('172')
+ elif band_n == 'B08' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('173')
+ elif band_n == 'B8A' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('174')
+ elif band_n == 'B09' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('175')
+ elif band_n == 'B10' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('176')
+ elif band_n == 'B11' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('177')
+ elif band_n == 'B12' and sensor.text == 'Sentinel-2A':
+ gscript.message(_(band_n))
+ text.write('178')
+ elif band_n == 'B01' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('179')
+ elif band_n == 'B02' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('180')
+ elif band_n == 'B03' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('181')
+ elif band_n == 'B04' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('182')
+ elif band_n == 'B05' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('183')
+ elif band_n == 'B06' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('184')
+ elif band_n == 'B07' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('185')
+ elif band_n == 'B08' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('186')
+ elif band_n == 'B8A' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('187')
+ elif band_n == 'B09' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('188')
+ elif band_n == 'B10' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('189')
+ elif band_n == 'B11' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('190')
+ elif band_n == 'B12' and sensor.text == 'Sentinel-2B':
+ gscript.message(_(band_n))
+ text.write('191')
+ else:
+ gscript.fatal('Bands do not seem to belong to a Sentinel image')
+ text.close()
+
+ if flags["a"]:
+ gscript.run_command('i.atcorr',
+ input=bb,
+ parameters=tmp_file,
+ output='{}_{}'.format(bb, suffix),
+ range='1,10000',
+ elevation=dem,
+ rescale=rescale,
+ flags='r')
+ cor_bands[key] = "{}_{}".format(bb, suffix)
+ else:
+ gscript.run_command('i.atcorr',
+ input=bb,
+ parameters=tmp_file,
+ output='{}_{}'.format(bb, suffix),
+ range='1,10000',
+ elevation=dem,
+ visibility=vis,
+ rescale=rescale,
+ flags='r')
+ cor_bands[key] = "{}_{}".format(bb, suffix)
+
+ gscript.message(_('--- All bands have been processed ---'))
+
+ if flags["t"]:
+ txt = open(txt_file, "w")
+ for key, value in cor_bands.items():
+ if str(key) in ['blue',
+ 'green',
+ 'red',
+ 'nir',
+ 'nir8a',
+ 'swir11',
+ 'swir12']:
+ txt.write(str(key) + '=' + str(value) + "\n")
+ txt.close()
+
+ for key, cb in cor_bands.items():
+ gscript.message(_(cb))
+ gscript.run_command('r.colors',
+ map=cb,
+ color='grey',
+ flags='e')
+
+ gscript.del_temp_region()
+ gscript.message(_('--- The computational region has been reset to the previous one ---'))
+
+if __name__ == "__main__":
+ options, flags = gscript.parser()
+ main()
+
Added: grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_ES.png
===================================================================
(Binary files differ)
Index: grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_ES.png
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_ES.png 2018-08-07 13:37:00 UTC (rev 73053)
+++ grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_ES.png 2018-08-07 13:42:24 UTC (rev 73054)
Property changes on: grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_ES.png
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property
Added: grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_GWF.png
===================================================================
(Binary files differ)
Index: grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_GWF.png
===================================================================
--- grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_GWF.png 2018-08-07 13:37:00 UTC (rev 73053)
+++ grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_GWF.png 2018-08-07 13:42:24 UTC (rev 73054)
Property changes on: grass-addons/grass7/imagery/i.sentinel.preproc/i_sentinel_preproc_GWF.png
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property
More information about the grass-commit
mailing list