[GRASS-SVN] r74080 - grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 12 02:30:03 PST 2019
Author: lucadelu
Date: 2019-02-12 02:30:03 -0800 (Tue, 12 Feb 2019)
New Revision: 74080
Modified:
grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.html
grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.py
Log:
i.sentinel.preproc: added topographic correction
Modified: grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.html
===================================================================
--- grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.html 2019-02-11 21:09:17 UTC (rev 74079)
+++ grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.html 2019-02-12 10:30:03 UTC (rev 74080)
@@ -1,12 +1,12 @@
<h2>DESCRIPTION</h2>
<em>i.sentinel.preproc</em> allows to import Sentinel-2
-images and perform atmospheric correction.
+images and perform atmospheric and topographic correction.
<p> <em>i.sentinel.preproc</em> 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>
+wraps the import, the atmospheric and the topographic correction using respectively
+<a href="i.sentinel.import.html">i.sentinel.import</a>,
+<a href="i.atcorr.html">i.atcorr</a> and <a href="i.topo.corr.html">i.topo.corr</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
@@ -33,6 +33,10 @@
corresponding code. In any case, <em>i.sentinel.preproc</em> writes a temporary control
file, changing it according to the band number, following the syntax rules and
codes of <em>i.atcorr</em> and then it runs <em>i.atcorr</em> for all bands.
+Using the <em>c</em> flag <em>i.sentinel.preproc</em> is able to perform also
+the topographic correction using <a href="i.topo.corr.html">i.topo.corr</a>
+creating the needed information as the illumination model based on the elevation
+model provided by the user.
<p>
<center>
@@ -272,6 +276,13 @@
into the corresponding code.
</ol>
+<h3>Topographic correction</h3>
+<em>i.sentinel.preproc</em> allows performing the topographic correction of all bands
+of a Sentinel-2 scene with a single process using <a href="i.topo.corr.html">i.topo.corr</a>.
+<em>i.sentinel.preproc</em> calculate the zenit and azimuth angles using
+<a href="r.sunmask.html">r.sunmask</a>, after that it create the illumination model based on
+the elevation model and apply it to all the bands of a Sentinel-2 scene
+
<h2>EXAMPLE</h2>
<p>
The example illustrates how to run <em>i.sentinel.preproc</em> for a Sentinel-2A image
Modified: grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.py
===================================================================
--- grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.py 2019-02-11 21:09:17 UTC (rev 74079)
+++ grass-addons/grass7/imagery/i.sentinel/i.sentinel.preproc/i.sentinel.preproc.py 2019-02-12 10:30:03 UTC (rev 74080)
@@ -92,6 +92,20 @@
#% required : no
#% guisection: Output
#%end
+#%option
+#% key: topo_method
+#% description: Topographic correction method
+#% options: cosine, minnaert, c-factor, percent
+#% required : no
+#% guisection: Input
+#%end
+#%option
+#% key: topo_prefix
+#% description: Prefix for topographic corrected images
+#% required : no
+#% answer: tcor
+#% guisection: Output
+#%end
#%flag
#% key: a
#% description: Use AOD instead visibility
@@ -112,6 +126,11 @@
#% description: Skip import of Sentinel bands
#% guisection: Input
#%end
+#%flag
+#% key: c
+#% description: Computes topographic correction of reflectance
+#% guisection: Input
+#%end
import grass.script as gscript
import xml.etree.ElementTree as et
@@ -130,7 +149,7 @@
bands = {}
cor_bands = {}
- dem = options['elevation']
+ dem = options['elevation']
vis = options['visibility']
input_dir = options['input_dir']
check_ndir = 0
@@ -140,7 +159,7 @@
level_dir = os.path.basename(input_dir).split('_')
# Check if the input directory is a .SAFE folder
if not input_dir.endswith('.SAFE'):
- gscript.fatal('The input directory is not a .SAFE folder. Please check the input directory')
+ gscript.fatal('The input directory is not a .SAFE folder. Please check the input directory')
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")]
@@ -165,7 +184,15 @@
processid = os.getpid()
txt_file = options['text_file']
tmp_file = gscript.tempfile()
+ topo_method = options['topo_method']
+ if topo_method and not flags['c']:
+ gscript.warning(_("To computes topographic correction of reflectance "
+ "please select also 'c' flag"))
+ elif flags['c'] and not topo_method:
+ gscript.warning(_("Topographic correction of reflectance will use "
+ "default method 'c-factor'"))
+
# Import bands
if not flags["i"]:
try:
@@ -182,7 +209,7 @@
# 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'):
@@ -470,7 +497,7 @@
text.write(str(25) + "\n")
elif sensor.text == 'Sentinel-2B':
text.write(str(26) + "\n")
- else:
+ else:
gscript.fatal('The input image does not seem to be a Sentinel image')
text.write('{} {} {:.2f} {:.3f} {:.3f}'.format(
time_py.month,
@@ -489,7 +516,7 @@
text.write('3' + "\n")
else: # Midlatitude summer
text.write('2' + "\n")
- elif lat < -15.00 and lat >= -45.00:
+ elif lat < -15.00 and lat >= -45.00:
if time_py.month in winter: # Midlatitude summer
text.write('2' + "\n")
else: # Midlatitude winter
@@ -642,7 +669,7 @@
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,
@@ -685,8 +712,43 @@
gscript.run_command('r.colors',
map=cb,
color='grey',
- flags='e')
+ flags='e', quiet=True)
+ if flags['c']:
+ gscript.message(_('--- Computes topographic correction of reflectance ---'))
+ dat = bb.split('_')[1]
+ # TODO understand better the timezone
+ sunmask = gscript.parse_command('r.sunmask', flags='sg', elevation=dem,
+ year=dat[0:4], month=int(dat[4:6]),
+ day=int(dat[6:8]), hour=int(dat[9:11]),
+ minute=int(dat[11:13]),
+ second=int(dat[13:15]), timezone=0)
+ z = 90. - float(sunmask['sunangleabovehorizon'])
+ if not topo_method:
+ topo_method = 'c-factor'
+ illu = "{}_{}_{}".format(bb, 'illu', processid)
+ gscript.run_command('i.topo.corr', flags='i', basemap=dem, zenit=z,
+ azimuth=sunmask['sunazimuth'], output=illu)
+ tcor = []
+ for ma in cor_bands.values():
+ out = "{}_double_{}".format(ma, processid)
+ tcor.append(out)
+ gscript.raster.mapcalc('{}=double({})'.format(out, ma))
+
+ gscript.run_command('i.topo.corr', basemap=illu, zenith=z,
+ input=','.join(tcor),
+ method=topo_method, output=options['topo_prefix'])
+ for ma in tcor:
+ inp = "{}.{}".format(options['topo_prefix'], ma)
+ gscript.run_command('g.rename', quiet=True,
+ raster="{},{}".format(inp,
+ inp.replace("_double_{}".format(processid),
+ "")))
+ gscript.run_command('g.remove', flags='f', type='raster', name=illu,
+ quiet=True)
+ gscript.run_command('g.remove', flags='f', type='raster',
+ name=','.join(tcor), quiet=True)
+
gscript.del_temp_region()
gscript.message(_('--- The computational region has been reset to the previous one ---'))
More information about the grass-commit
mailing list