GRASS-Wiki - User contributions [en]

GRASS Community Sprint Bonn 2016

2016-02-25T00:42:18Z

⚠️Andbal: corrected host country

{{toc|right}}
During the [http://2016.foss4g.org/ FOSS4G 2016] , 21-28 August 2016, [http://en.wikipedia.org/wiki/Bonn Bonn], [https://en.wikipedia.org/wiki/Germany Germany], the GRASS GIS project will join the code sprint

'''Important dates:'''

* CodeSprint I: Sun 21.08.2016 12:00 - Tue 23.08.2016 (open end)
* [http://2016.foss4g.org/#page9 FOSS4G 2016 Schedule]
* CodeSprint II: Sat 27.08.2016 - Sun 28.08.2016 12:00 a.m

== Purpose ==

Get together of GRASS GIS project members to take decisions and tackle small and larger problems.

For this meeting, we welcome people committed to improving the GRASS GIS project. This includes developers, documenters, bug reporters, translators, and users!

== Sponsors ==

If you are interested to sponsor the GRASS Community Sprint, please read about 
    '''sponsoring the GRASS project at [http://grass.osgeo.org/donations/ http://grass.osgeo.org/donations/]''' 
and, if needed, contact [[User:Neteler|Markus Neteler]] <tt><neteler at osgeo.org></tt>. Any surplus at the end of the event will be turned over to the GRASS GIS project.

The GRASS Community Sprint is a great occasion for you to support the development of GRASS. With your contribution you'll enable more developers to meet in Prague. The community sprint is an important opportunity for the GRASS developers to discuss and collaboratively resolve bugs, plan the direction for the project and work on new features. Please see below for the more detailed agenda. The developers and contributors are donating their valuable time, so it would be great if in-kind funding can be made available from within the community to cover out-of-pocket expenses. All of the work that takes place at the community sprint will be directly contributed back into the GRASS project to the benefit of everyone who uses it.

=== Thanks to our sponsors ===

== Venue ==

* Bonn, Germany: http://www.basecamp-bonn.de

== Accommodation and Costs ==

* to be defined

== Weather ==

* https://en.wikipedia.org/wiki/Bonn#Climate

== Agenda - What we plan to do ==

'''Note:''' The program is generally open for your ideas. Please edit this wiki page!

== Participation ==

For organizationals reasons, please '''additionally register''' at: https://wiki.osgeo.org/wiki/FOSS4G_2016_Code_Sprint

=== In person ===

Please add your name here or contact [[User:Lucadelu|Luca Delucchi]] <tt><lucadeluge at gmail com></tt>:

{|class="wikitable sortable" border="2" cellspacing="0" cellpadding="4" rules="all" style="margin:1em 1em 1em 0; border:solid 1px #AAAAAA; border-collapse:collapse; background-color:#edf9c7; font-size:95%; empty-cells:show;"
|width=50px|'''No'''
|width=130px|'''Participant '''
|width=100px|'''Country'''
|width=100px|'''Arrival'''
|width=100px|'''Departure'''
|'''Topic'''
|width=75px|'''T-Shirt'''
|'''Notes'''
|-
|1
|[[User:Lucadelu|Luca Delucchi]]
|Italy
| August 21
| August 28
|
| L
|
|-
|2
|[[User:Landa|Martin Landa]]
|Czech Republic
| August 21
| August 28
|
| L
|
|-
|3
|[[User:Neteler|Markus Neteler]]
|Italy
| August 21
| August 28
|
| M
|
|-
|4
|[[User:Ychemin|Yann Chemin]]
|France
| August 21
| August 28
|
| XL
|
|-
|5
|[[User:NikosA|Nikos Alexandris]]
|Switzerland
| August 21
| August 28
|
| S or M
|
|-
|}

=== Via IRC chat ===

== Individual Preparation ==

* Bring your own computer
* Bring [http://en.wikipedia.org/wiki/Mains_electricity_by_country your power connector adapter] if needed ([https://en.wikipedia.org/wiki/File:Schuko_plug_and_socket.png Germany])
* Install subversion and the compiler tools, and come with a working GRASS development environment if possible.

== Photos ==

We have to remember to take at least one group photo! :-)

== FAQ ==

* ''Is the GRASS Community Sprint just a coding event?''
** It is mainly a coding and documentation event. It is a working session for people who are already participants in the GRASS project and/or are committed to improving the GRASS project.
** On demand we can do some presentations of current working GRASS implementation and new upcoming features to spread the idea of Open Source GIS software
* ''Is the GRASS Community Sprint for developers only?''
** No: anybody can help, with testing, checking out bugs and fixes, documentation and more.
* ''Where can I get help and more information about the community sprint?''
** Contact [[User:Lucadelu|Luca Delucchi]] <tt><lucadeluge at gmail com></tt>

[[Category: Workshops]]
[[Category: Code Sprint]]
[[Category: 2016]]

LANDSAT

2016-02-16T22:19:31Z

⚠️Andbal: added info on haze removal

== Data Availability==

Landsat imagery can be obtained from [http://landsat.usgs.gov/ USGS' LANDSAT archive]:

* [http://glovis.usgs.gov USGS GloVis]
* [http://earthexplorer.usgs.gov/ EarthExplorer]
* [http://earthexplorer.usgs.gov/bulk/help Bulk Download Orders] using the [http://earthexplorer.usgs.gov/bulk/ Bulk Download Application]
* [http://landsatlook.usgs.gov/ LandsatLook Viewer]

See also the [[Global datasets]] wiki page

On 06/12/13, USGS [http://landsat.usgs.gov/about_LU_Vol_7_Issue_4.php announced] the availability of [http://landsat.usgs.gov/landsat8.php LANDSAT 8]. See also [http://landsat.usgs.gov/band_designations_landsat_satellites.php this page].

== Modules overview ==

GRASS GIS features a complete set of native modules and various add-ons for pre- and post-processing Landsat satellite imagery. The following lists offer an overview of related modules and add-ons.

=== Generic modules applicable to Landsat ===

* {{cmd|d.rgb}} → display 3-band data
* {{cmd|r.composite}} → flatten 3-bands of data into a single image (lossy, maybe used in combination with {{cmd|i.landsat.rgb|version=64}} (GRASS 6.x) | {{cmd|i.colors.enhance|version=70}} (GRASS 7.x)
* {{cmd|i.atcorr}} → correct top of atmosphere to surface reflectance - see also the [[Atmospheric correction]] wiki page
* {{cmd|i.topo.corr}} → used to topographically correct reflectance from imagery files, e.g. obtained with {{cmd|i.landsat.toar}}, using a sun illumination terrain model
* {{cmd|i.oif}} → calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
* {{AddonCmd|i.histo.match}} (addon) → histogram matching of two or more raster maps (in grass 7) ''Note, the module works with integer values and does not accept the "." character as part of the raster map's name!''

=== Landsat specific modules ===

* {{cmd|i.landsat.rgb|version=64}} (GRASS 6.x) | {{cmd|i.colors.enhance|version=70}} (GRASS 7.x) → auto-enhance colors
* {{cmd|i.landsat.toar}} → convert DN to top of atmosphere radiance
* {{cmd|i.landsat.acca}} → cloud cover assessment

=== Landsat specific GRASS AddOns ===

* {{AddonCmd|i.landsat.trim|version=64}} → trims border fringes for each band separately or with the MASK where coverage exists for all bands
* {{AddonCmd|i.landsat.dehaze|version=64}} (addon) → haze removal

== Pre-Processing ==

=== Overview ===

Typically, pre-processing Landsat imagery comprises the following steps:

# '''import''' in the database → {{cmd|r.in.gdal}}
# geometrically & orthometrically correct imagery
#* already done for L1T products, read more at [http://landsat.usgs.gov//Landsat_Processing_Details.php USGS' Landsat Information Products] webpage
# optionally, automatically '''cut-off border fringes''' → {{AddonCmd|i.landsat.trim|version=64}}
#* of course one can use the official WRS2 Path/Row vector tiles to manually trim border fringes
# optionally, '''denoise''' for obvious/intensive salt & pepper effects, stripes, etc.
#* for example by applying Principal Components Analysis as a denoising technique → {{cmd|i.pca}} ''<<< Re-order this step?''
# '''convert''' the '''Digital Numbers''' (DNs) to '''Top-of-Atmosphere Radiances/Reflectances''' (ToARs) → {{cmd|i.landsat.toar}}
# optionally, '''correct for atmospheric effects''' → {{cmd|i.atcorr}}
#* that is, accounting for distorting atmospheric effects and estimating actual reflectances as they would have been measured on the ground
#* also described as conversion to Top-of-Canopy Reflectances (ToCRs)
# '''assessing cloud cover''' → {{cmd|i.landsat.acca}}
#* optionally, detect and remove clouds shadows as well
# '''topographically normalise''' imagery → {{cmd|i.topo.corr}}
#* also known as topographic correction, that is, accounting for illumination differences due to the acquisition's geometry
# '''radiometrically normalise''' → one approach via {{AddonCmd|i.histo.match}} (in '''grass 7''', addon), also known as relative radiometric normalisation -- one approach is the ''histogram matching'' technique of two or more raster maps
#* '''ToDo/More techniques?'''

See also
* [http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ Processing Landsat 8 data in GRASS GIS 7: Import and visualization]

=== Importing data ===

Importing Landsat spectral bands in GRASS GIS' data base can be done both from the Graphical User Interface or from the command line

# Open GRASS GIS, select ''Location Wizard'' in order to create a new location from georeferenced file
# Use the ''Import file tool'', or the {{cmd|r.in.gdal}} module, to import the GeoTIFF files into GRASS GIS.

See also [http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import Automated data import] below.

==== Notes ====

# Most Landsat scenes are delivered in ''north-is-up'' orientation, hence the import process is straightforward.
# If you get <source lang="bash" enclose="none">ERROR: Input map is rotated - cannot import.</source>, the image must be first rotated externally, applying the rotation info stored in the metadata field of the raster image file. For example, the <code>gdalwarp</code> software can be used to transform the map to North-up (note, there are several <code>gdalwarp</code> parameters to select the resampling algorithm):

<ul>
<source lang="bash">
gdalwarp rotated.tif northup.tif
</source>
</ul>

==== Hint: Minimal disk space copies ====

Here's a little trick using {{cmd|r.reclass}} to rename maps (for example, from ''L71074092_09220040924_B'''10''''' to ''L71074092_09220040924_B'''1''''') without touching the data or wasting disk space:

<source lang="bash">
BASE=L71074092_09220040924

for BAND in 10 20 30 40 50 61 70 80; do
BAND1st=`echo $BAND | sed -e 's/0$//'`
r.reclass in="${BASE}_B$BAND" out=$BASE.$BAND1st << EOF
* = *
EOF
done
</source>

==== Automated data import ====

The following ''Python'' script imports Landsat imagery into GRASS' data base. Specifically, the script

* creates an independent Mapset for each Landsat scene
* imports and renames bands of a scene as <tt>B<id></tt>, e.g. B10, B20, ..., B80.
* additionaly sets up the timestamp based on MTL metadata file

'''Note,''' the (newest) official naming pattern for Landsat scenes -- explained in [https://lta.cr.usgs.gov/landsat_dictionary.html USGS' Landsat Data Dictionary] as the [https://lta.cr.usgs.gov/landsat_dictionary.html#entity_id Landsat Scene Identifier] -- and all individual bands that compose a scene -- bands have a suffix which is like _B10, _B20, _B30, etc. -- differ from what some Landsat specific modules expect. For example, the modules {{cmd|i.landsat.toar}} and {{cmd|i.landsat.acca}} expect the bands to follow a naming pattern such as "scenename.1, .2, .3", etc.

To use the script save it as <source lang="bash" enclose="none">import_landsat.py</source> file and make sure it is granted the execution permission.

<source lang=python>
#!/usr/bin/python

import os
import sys
import glob
import grass.script as grass

def get_timestamp(mapset):
try:
metafile = glob.glob(mapset + '/*MTL.txt')[0]
except IndexError:
return

result = dict()
try:
fd = open(metafile)
for line in fd.readlines():
line = line.rstrip('\n')
if len(line) == 0:
continue
if any(x in line for x in ('DATE_ACQUIRED', 'ACQUISITION_DATE')):
result['date'] = line.strip().split('=')[1].strip()
finally:
fd.close()

return result

def import_tifs(mapset):
meta = get_timestamp(mapset)
for file in os.listdir(mapset):
if os.path.splitext(file)[-1] != '.TIF':
continue
ffile = os.path.join(mapset, file)
if ('VCID') in ffile:
name = "".join((os.path.splitext(file)[0].split('_'))[1::2])
else:
name = os.path.splitext(file)[0].split('_')[-1]
if len(name) == 3 and name[-1] == '0':
band = int(name[1:2])
elif len(name) == 3 and name[-1] != '0':
band = int(name[1:3])
else:
band = int(name[-1:])
grass.message('Importing %s -> %s@%s...' % (file, name, mapset))
grass.run_command('g.mapset',
flags = 'c',
mapset = mapset,
quiet = True,
stderr = open(os.devnull, 'w'))
grass.run_command('r.in.gdal',
input = ffile,
output = name,
quiet = True,
overwrite = True,
title = 'band %d' % band)
if meta:
year, month, day = meta['date'].split('-')
if month == '01':
month = 'jan'
elif month == '02':
month = 'feb'
elif month == '03':
month = 'mar'
elif month == '04':
month = 'apr'
elif month == '05':
month = 'may'
elif month == '06':
month = 'jun'
elif month == '07':
month = 'jul'
elif month == '08':
month = 'aug'
elif month == '09':
month = 'sep'
elif month == '10':
month = 'oct'
elif month == '11':
month = 'nov'
elif month == '12':
month = 'dec'

grass.run_command('r.timestamp',
map = name,
date = ' '.join((day, month, year)))

def main():
if len(sys.argv) == 1:
for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
import_tifs(directory)
else:
import_tifs(sys.argv[1])

if __name__ == "__main__":
main()
</source>

'''Example of usage'''

After having collected the Landsat scenes of interest in one directory, the script can be used as follows:

* <source lang="python" enclose="none">./import_landsat.py</source> → the script will walk through a ''pool'' directory that contains unique Landsat scene directories (e.g. three directories named after the official naming pattern: LT51800342011158MOR00 LT51810352009079MTI00 LT51820352009326MTI00) and import all bands of each individual scene in their own Mapset
* <source lang="python" enclose="none">./import_landsat.py LM41890261983200FFF03</source> → the scrip will import only bands of the specified Landsat scene directory

== Post-Processing ==

=== Natural color composites ===

Creating natural (also known as true-) color composites, can be done by

a) using {{cmd|i.landsat.rgb|version=64}} (GRASS 6.x) | {{cmd|i.colors.enhance|version=70}} (GRASS 7.x) to automatically balance the colors for the Red, Green and Blue bands

or

b) equalizing the colors on each of the Red, Green and Blue bands

<ul>
<source lang="bash">
r.colors -e map=band1 color=grey
</source>
</ul>

'''Composites''' can then be produced with the {{cmd|r.composite}} module.

'''Resetting the color tables''' of all bands back to normal greyscale can be done with a for loop:

<source lang="bash">
BASE=L71074092_09220040924

for map in `g.mlist pat="$BASE.[0-8]*"` ; do
r.colors $map color=grey255
done
</source>

See also:
* [http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening]

=== Create a MASK to only show data where coverage exists for all bands ===

<source lang="bash">
BASE=L71074092_09220040924

g.region rast=$BASE.1
r.series in=`g.mlist pat="$BASE.[0-8]*" sep=,` -n out=$BASE.thresh method=threshold thresh=1
r.mapcalc "$BASE.mask = if(isnull($BASE.thresh))"
g.remove $BASE.thresh
</source>

=== Calculate Top-of-Atmosphere Reflectance and band-6 Temperature ===

Calculate ''Top-of-Atmosphere'' reflectance and band-6 temperature using the {{cmd|i.landsat.toar}} module. For details, refer to USGS' [http://landsat.usgs.gov/Landsat_Metadata_Changes.php Landsat Filename and Metadata Changes] dedicated page. For Landsat 8, see also [http://landsat.usgs.gov/about_LU_Vol_7_Issue_4.php#3a here].
<source lang="bash">
i.landsat.toar input_prefix=$BASE output_prefix=${BASE}_toar metfile=${BASE}_MTL.txt
</source>

Note, the resulting temperature map is in Kelvin:
<source lang="bash">
# convert to degree Celsius
r.mapcalc "$BASE.temp_celsius = ${BASE}_toar.6 - 273.15"
r.info -r $BASE.temp_celsius
</source>

=== Haze removal ===

Simple haze removal can be done with {{AddonCmd|i.landsat.dehaze|version=64}}. This addons applies a bandwise haze correction using tasscap4 (haze) and linear regression.

During the 2000s, prior acquired Landsat data were reprocessed to LPGS (Level 1 Product Generation System). Seems that with this level of processing, haze effect in raw data was removed and now is sufficient to apply {{Cmd|i.colors.enhance}} to obtain a good image with high contrast.

=== Atmospheric correction ===

See [[Atmospheric correction]]

=== Cloud identification ===

Identify clouds in the image with {{cmd|i.landsat.acca}}:
<source lang="bash">i.landsat.acca -f input_prefix=226_62_toar. output=226_62.acca_cloudmask</source>

Mask out the clouds:
<source lang="bash">r.mapcalc "MASK = if(isnull($BASE.acca_cloudmask))"</source>

== Download sample data ==

=== Preprocessed Landsat-7 data for North Carolina ===

The [[Sample_datasets|North Carolina 2008 sample dataset]] comes with 3 different Landsat scenes:
: (''Wake County -- path: 16 row: 35 for various dates'')

The above import efforts are not needed since the data are already in a GRASS location.

=== Landsat-5: Oct 14, 1987 ===

;Glovis download: LT50160351987287XXX08
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LT50160351987287XXX08&dataset_name=LANDSAT_TM



Provided metadata:
IMAGE_ID=P016R35_5T871014
PATH=16
ROW=35
DATE=10/14/87
PLATFORM=LANDSAT5

=== Landsat-7: Mar 31, 2000 ===

;Glovis download: LE70160352000091EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352000091EDC00&dataset_name=LANDSAT_ETM

Values in the metadata below indicate that the version provided with the NC 2008 dataset's production date was after July 1, 2000, and that the channel gains were <tt>HHHLHLHHL</tt>.

Convert DNs to radiance/temperatures:<source lang="bash">
GRASS> i.landsat.toar -v band=lsat7_2000 sensor=7 date=2000-03-31 \
product_date=2000-07-02 solar_elevation=51.5246529 gain=HHHLHLHHL</source>

Provided metadata:
SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
'''ACQUISITION_DATE=2000-03-31'''
WRS_PATH=16
CPF_FILE_NAME=L7CPF20000101_20000331_12
SUN_AZIMUTH=139.6033279
'''SUN_ELEVATION=51.5246529'''
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

=== Landsat-7: May 24, 2002 ===

;Glovis download: LE70160352002144EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352002144EDC00&dataset_name=LANDSAT_ETM



Provided metadata: (`p016r035_7x20020524.met`)
SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
'''ACQUISITION_DATE=2002-05-24'''
WRS_PATH=016
WRS_ROW=035
SUN_AZIMUTH=120.8810347
'''SUN_ELEVATION=64.7730999'''
QA_PERCENT_MISSING_DATA=0
CLOUD_COVER=0
CPF_FILE_NAME=L7CPF20020401_20020630_03
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

== LANDSAT Image classification ==

See [[Image classification]]

== Time series analysis ==

See [[Time series]]

== Sources ==

* [http://landsat.usgs.gov/tools_access_all_faqs.php Frequently Asked Questions about the Landsat Missions]
* [https://landsat8portal.eo.esa.int/portal/ European LANDSAT 8 data in near-realtime]
* [http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ Processing Landsat 8 data in GRASS GIS 7: Import and visualization]
* [http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening]

[[Category: Documentation]]
[[Category: FAQ]]
[[Category: Landsat]]
[[Category: Image processing]]
[[Category: Import]]
[[Category: Digital Numbers]]
[[Category: Radiance]]
[[Category: Reflectance]]
[[Category: Geodata]]

GRASS and Python

2016-01-26T10:33:49Z

⚠️Andbal: update link to trac wiki

''This page discusses usage GRASS GIS and Python in a general way if you just want to write a script to run in GRASS GIS go to [[GRASS Python Scripting Library]] wiki page or [http://grass.osgeo.org/grass70/manuals/libpython/ Python API for GRASS GIS 7] documentation.''

==Writing Python scripts in GRASS==

Python is a programming language which is more powerful than shell scripting but easier and more forgiving than C.
The Python script can contain simple module description definitions which will be processed with {{cmd|g.parser}}, as shown in the example below. In this way with no extra coding a GUI can be built, inputs checked, and a skeleton help page can be generated automatically. In addition it adds links to the GRASS message translation system. The library for "scripting" is "grass.script", typically used as:

<source lang="python">
import grass.script as gscript
</source>

The related files are at $GISBASE/etc/python/grass/script/*.py. See below for more details.

''Note: For code which needs access to the power of C, you can access the GRASS C library functions via [[GRASS and Python#Python Ctypes Interface|the Python "ctypes" interface]].''

=== Integrated Python script editor ===

The [[wxGUI]] Layer Manager in GRASS 6.4.3+ comes with a "Python shell" which enables users to type and execute python commands directly in wxGUI environment.

[[Image:wxgui-pyshell.png|center|400px|Embedded interactive Python Shell in wxGUI Layer Manager]]

=== External Python editors (IDE) ===

Besides the wxGUI Python script editor also advanced editors can be used like '''Spyder''' (The Scientific PYthon Development EnviRonment). For details, see [[Tools for Python programming]]

=== Using the GRASS Python Scripting Library ===

You can run Python scripts easily in a GRASS session.

To write these scripts,
* check the code in lib/python/ which provides grass.script in order to support GRASS scripts written in Python. 
See the [[GRASS Python Scripting Library]] for notes and examples.
* The scripts/ directory of GRASS contains a series of examples actually provided to the end users.

For the desired Python code style, have a look at [https://trac.osgeo.org/grass/wiki/Submitting/Python Submitting Python].

=== Creating Python scripts that call GRASS functionality from outside ===

Note: This is a more advanced use case of using GRASS' functionality from outside via Python. Commonly, a user will run GRASS Python script from inside a GRASS session, i.e. either from the command line or from the Python shell embedded in the wxGUI ([[:File:Wxgui-pyshell.png|screenshot]]).

For calling GRASS functionality from outside, see also [[Working with GRASS without starting it explicitly]].

The next two subsections (MS-Windows and Linux) are somewhat outdated, to be improved (see grass-dev archives, 5/2014

==== MS-Windows ====

In order to use GRASS functionality via Python from outside, some environment variables have to be set:

GISBASE= C:\GRASS-64
GISRC= C:\Documents and Settings\user\.grassrc6
LD_LIBRARY_PATH= C:\GRASS-64\lib
PATH= C:\GRASS-64\etc;C:\GRASS-64\etc\python;C:\GRASS-64\lib;C:\GRASS-64\bin;C:\GRASS-64\extralib;C:\GRASS-64\msys\bin;C:\Python26;
PYTHONLIB= C:\Python26
PYTHONPATH= C:\GRASS-64\etc\python
GRASS_SH= C:\GRASS-64\msys\bin\sh.exe

Some hints:

# The ".grassrc6" file listed above must exist. Run GRASS one time interactively to create it or see below for how to writing it
# The Python interpreter (python.exe) needs to be in the PATH
# Python needs to be associated with the .py extension
# PATHEXT needs to include .py if you want to be able to omit the extension
# PYTHONPATH needs to be set to %GISBASE%\etc\python

Points 2-4 should be taken care of by the Python installer. 5 needs to be done by the startup (currently, this doesn't appear to be the case on MS-Windows).

Alternatively to run GRASS interactively, you can also create the ".grassrc6" file yourself, e.g. (update to existing directory for "grassdata"):

GISDBASE: C:\Documents and Settings\user\grassdata
LOCATION_NAME: nc_spm_08
MAPSET: user1
GRASS_DB_ENCODING: ascii

It doesn't matter what the file is called, so long as %GISRC% points to it and it contains the necessary settings.

The normal location for GRASS 6.x on Windows is:

%APPDATA%\GRASS6\grassrc6

On Windows 7, a typical setting for %APPDATA% is

C:\Users\<username>\AppData\Roaming

==== Linux ====

In order to use GRASS functionality via Python from outside, some environment variables have to be set:

export GISBASE="/usr/local/grass-6.4.svn/"
export PATH="$PATH:$GISBASE/bin:$GISBASE/scripts"
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$GISBASE/lib"
# for parallel session management, we use process ID (PID) as lock file number:
export GIS_LOCK=$$
# path to GRASS settings file
export GISRC="$HOME/.grassrc6"
export PYTHONPATH="$PYTHONPATH:$GISBASE/etc/python"

=== Executing GRASS Python scripts without starting GRASS explicitly ===

* See [[Working_with_GRASS_without_starting_it_explicitly#Python:_GRASS_GIS_7|Working with GRASS without starting it explicitly]] (examples for GRASS GIS 6 and 7)

=== Running external commands from Python ===
For information on running external commands from Python, see:
http://docs.python.org/lib/module-subprocess.html

Avoid using the older os.* functions. Section 17.1.3 lists equivalents
using the Popen() interface, which is more robust (particularly on
Windows).

=== Testing and installing Python extensions ===

==== Debugging ====

Make sure the script is executable:

chmod +x /path/to/my.extension.py

During development, a Python script can be debugged using the Python Debugger (pdb):

python -m pdb /path/to/my.extension.py input=my_input_layer output=my_output_layer option=value -f

==== Installation ====

Once you're happy with your script, you can put it in the scripts/ folder of your GRASS install. To do so, first create a directory named after your extension, then create a Makefile for it, and a HTML man page:

cd /path/to/grass_src/
cd scripts
ls # It is useful to check out the existing scripts and their structure
mkdir my.extension
cd my.extension
cp path/to/my.extension.py .
touch my.extension.html
touch Makefile

Next step is to edit the Makefile. It is a very simple text file, the only thing to check is to put the right extension name (WITHOUT the .py file extension) after PGM:

MODULE_TOPDIR = ../..

PGM = my.extension

include $(MODULE_TOPDIR)/include/Make/Script.make

default: script

The HTML file would be generated automatically. If you want to add more precisions in it, you can do it (just make sure you start at DESCRIPTION. See existing scripts.)

You can then run "make" within the my.extension folder. Running "make" in the extension directory places the resulting files in the staging directory (path/to/grass_src/dist.<YOUR_ARCH>/). If you're running GRASS from the staging directory (/path/to/grass_src/bin.<YOUR_ARCH>/grass7), subsequent commands will used the updated files.

# in your extension directory (/path/to/grass_src/scripts/my.extension/)
make
# Starting GRASS from the staging directory
/path/to/grass_src/bin.<YOUR_ARCH>/grass7
my.extension help

You can also run "make install" from the top level directory of your GRASS install (say /usr/local/src/grass_trunk/). Running "make install" from the top level just copies the whole of the dist.<YOUR_ARCH>/ directory to the installation directory (e.g. /usr/local/grass70) and the bin.<YOUR_ARCH>/grass70 bin file to the bin directory (e.g. /usr/local/bin), and fixes any embedded paths in scripts and configuration files.

cd /path/to/grass_src
make install
# Starting GRASS as usual would work and show your extension available
grass7
my.extension help

=== Using the parser ===

Be sure to use the standardized option G_OPT_*_INPUT (e.g. G_OPT_R_INPUT), the definitions can be found here:
* http://grass.osgeo.org/programming7/parser__standard__options_8c_source.html
* Overview table of [https://grass.osgeo.org/grass71/manuals/parser_standard_options.html parser standard options]

For more details, see
* [[Module command line parser]]
* [[GRASS_Python_Scripting_Library#Parsing_the_options_and_flags]]

=== Using a version of Python different from the default installation ===

Sometimes you don't want to use the default Python version installed on your system. You can control the interpreter used for building by overriding the PYTHON make variable, e.g.

make PYTHON=/path/to/python ...

For controlling the interpreter used at run time:
* On Windows, Python scripts are invoked via %GRASS_PYTHON%, so changing that environment variable will change the interpreter.
* On Unix, GRASS_PYTHON is only used for components which use wxPython, e.g. wxGUI or the parameter dialogs which are displayed when modules are run without arguments.

Scripts rely upon the "#!/usr/bin/env python" line. You can control
the interpreter used by scripts by creating a directory containing a
symlink named "python" which refers to your preferred version of the
interpreter, and placing that directory at the front of $PATH.

==Python extensions in GRASS GIS==

=== Python Scripting Library ===

* See [[GRASS Python Scripting Library]]

=== pygrass Library ===

Pygrass is an object oriented Python Application Programming Interface (API) for GRASS GIS. pygrass wants to improve the integration between GRASS and Python, render the use of Python under GRASS more consistent with the language itself and make the GRASS scripting and programming activity easier and more natural to the final users. Some simple examples below.

For more details, see '''[[Python/pygrass|pygrass]]'''

'''Example 1''':
<source lang="python">
from grass.pygrass.modules import Module

slope_aspect = Module("r.slope.aspect")
slope_aspect(elevation='elevation', slope='slp', aspect='asp',
format='degrees', overwrite=True)
</source>

'''Example 2''', using "shortcuts":
<source lang="python">
#!/usr/bin/env python

# simple example for pyGRASS usage: raster processing via modules approach

from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r

g.message("Filter elevation map by a threshold...")

# set computational region
input = 'elevation'
g.region(rast=input)

# hardcoded:
# r.mapcalc('elev_100m = if(elevation > 100, elevation, null())', overwrite = True)
# with variables

output = 'elev_100m'
thresh = 100.0
r.mapcalc("%s = if(%s > %d, %s, null())" % (output, input, thresh, input), overwrite = True)
r.colors(map=output, color="elevation")
</source>

'''Example 3''', using "shortcuts" and external raster format support:
<source lang="python">
#!/usr/bin/env python

# simple example for pyGRASS usage: raster processing via Modules approach
# Read GeoTIFF directly, write out GeoTIFF directly
import os
import tempfile
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules import Module

# use alias name for GRASS GIS commands with more than one dot in name
r.external_out = Module('r.external.out')

# get the raster elevation map into GRASS without import
geotiff = 'elevation.tif' # with path
curr_raster='myraster'
r.external(input=geotiff, output=curr_raster)

# define output directory for files resulting from GRASS calculation:
gisdb = os.path.join(tempfile.gettempdir(), 'gisoutput')
try:
os.stat(gisdb)
except:
os.mkdir(gisdb)
g.message("Output will be written to %s" % gisdb)

# define where to store the resulting files and format
r.external_out(directory=gisdb, format='GTiff')

# do the GRASS GIS job
g.message("Filter elevation map by a threshold...")
# set computational region
g.region(rast=curr_raster)

# generate GeoTIFF
output = 'elev_100m.tif'
thresh = 100.0
r.mapcalc("%s = if(%s > %d, %s, null())" % (output, curr_raster, thresh, curr_raster), overwrite = True)
r.colors(map=output, color="elevation")

# cease GDAL mode, back to GRASS data writing mode
r.external_out(flags = 'r')
# the resulting GeoTIFF file is in the gisdb directory
</source>

'''Example 4''':
<source lang="python">
#!/usr/bin/env python

# Example for pyGRASS usage - vector API

from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.vector import VectorTopo
from grass.pygrass.modules import Module

vectmap = 'myzipcodes_wake'

g.message("Importing SHAPE file ...")
ogrimport = Module('v.in.ogr')
ogrimport('/home/neteler/gis_data/zipcodes_wake.shp', output=vectmap)

g.message("Assessing vector topology...")
zipcodes = VectorTopo(vectmap)

# Open the map with topology:
zipcodes.open()

# query number of topological features
areas = zipcodes.number_of("areas")
islands = zipcodes.number_of("islands")
print 'Map: <' + vectmap + '> with %d areas and %d islands' % (areas, islands)

dblink = zipcodes.dblinks[0]
print 'DB name:'
print dblink.database
table = dblink.table()
print 'Column names:'
print table.columns.names()
print 'Column types:'
print table.columns.types()

zipcodes.close()
</source>

* For more examples, see '''[[Python/pygrass|pygrass]]'''

=== Python Ctypes Interface ===

This interface allows calling GRASS library functions from Python scripts. See [[Python Ctypes Examples]] for details.

Examples:

* GRASS 7: [http://trac.osgeo.org/grass/browser/grass/trunk/doc/python/raster_example_ctypes.py raster], [http://trac.osgeo.org/grass/browser/grass/trunk/doc/python/vector_example_ctypes.py vector] example

* Latest and greatest: GRASS 7 Python {{src|scripts}}

* <strike>[[PythonSwigExamples|More complicated examples]]</strike> <<-- TODO: update to Ctypes

Sample script for GRASS 6 raster access (use within GRASS, Spearfish session):
<source lang="python">
#!/usr/bin/env python

## TODO: update example to Ctypes

import os, sys
from grass.lib import grass

if "GISBASE" not in os.environ:
print "You must be in GRASS GIS to run this program."
sys.exit(1)

if len(sys.argv)==2:
input = sys.argv[1]
else:
input = raw_input("Raster Map Name? ")

# initialize
grass.G_gisinit('')

# find map in search path
mapset = grass.G_find_cell2(input, '')

# determine the inputmap type (CELL/FCELL/DCELL) */
data_type = grass.G_raster_map_type(input, mapset)

infd = grass.G_open_cell_old(input, mapset)
inrast = grass.G_allocate_raster_buf(data_type)

rown = 0
while True:
myrow = grass.G_get_raster_row(infd, inrast, rown, data_type)
print rown, myrow[0:10]
rown += 1
if rown == 476:
break

grass.G_close_cell(inrast)
grass.G_free(cell)
</source>

Sample script for vector access (use within GRASS, Spearfish session):

<source lang="python">
#!/usr/bin/python

# run within GRASS Spearfish session
# run this before starting python to append module search path:
# export PYTHONPATH=/usr/src/grass70/swig/python
# check with "import sys; sys.path"
# or:
# sys.path.append("/usr/src/grass70/swig/python")
# FIXME: install the grass bindings in $GISBASE/lib/ ?

import os, sys
from grass.lib import grass
from grass.lib import vector as grassvect

if "GISBASE" not in os.environ:
print "You must be in GRASS GIS to run this program."
sys.exit(1)

if len(sys.argv)==2:
input = sys.argv[1]
else:
input = raw_input("Vector Map Name? ")

# initialize
grass.G_gisinit('')

# find map in search path
mapset = grass.G_find_vector2(input,'')

# define map structure
map = grassvect.Map_info()

# define open level (level 2: topology)
grassvect.Vect_set_open_level (2)

# open existing map
grassvect.Vect_open_old(map, input, mapset)

# query
print 'Vect map: ', input
print 'Vect is 3D: ', grassvect.Vect_is_3d (map)
print 'Vect DB links: ', grassvect.Vect_get_num_dblinks(map)
print 'Map Scale: 1:', grassvect.Vect_get_scale(map)
print 'Number of areas:', grassvect.Vect_get_num_areas(map)

# close map
grassvect.Vect_close(map)
</source>

=== Using ipython ===

Using [http://ipython.org/ ipython] is great due to the TAB completion. Don't remember the function name or curious about what's offered? Hit TAB (maybe twice for a list).

Here an example how to connect to C functions directly from Python:

<source lang="python">
grass70
# ... in the GRASS GIS terminal
ipython
In [1]: from grass.lib import gis

In [2]: gis.<TAB>
Display all 729 possibilities? (y or n)
gis.ARRAY gis.G_fatal_longjmp gis.G_set_key_value
gis.ArgumentError gis.G_file_name gis.G_set_keywords
gis.Array gis.G_file_name_misc gis.G_set_ls_exclude_filter
...
gis.G_fatal_error gis.G_set_gisrc_mode gis.wstring_at

In [2]: gis.G_message("Hello world")
Hello world
</source>

A full ipython workshop is available here: https://github.com/wenzeslaus/python-grass-addon

=== wxPython GUI development ===

* See the [[wxGUI]] wiki page

=== Python-GRASS add-ons ===

Stand-alone addons:

* Jáchym Čepický's PyWPS, GRASS-Web Processing Service (http://pywps.wald.intevation.org)

=== Using GRASS gui.tcl in Python ===

Here is some example code to use the grass automatically generated guis in python code. This could (should) all be bundled up and abstracted away so that the implementation can be replaced later.

<source lang="python">
import Tkinter
import os

# Startup (once):

tk = Tkinter.Tk()
tk.eval ("wm withdraw .")
tk.eval ("source $env(GISBASE)/etc/gui.tcl")
# Here you could do various things to change what the gui does
# See gui.tcl and README.GUI

# Make a gui (per dialog)
# This sets up a window for the command.
# This can be different to integrate with tkinter:
tk.eval ('set path ".dialog$dlg"')
tk.eval ('toplevel .dialog$dlg')
# Load the code for this command:
fd = os.popen ("d.vect --tcltk")
gui = fd.read()
# Run it
tk.eval(gui)
dlg = tk.eval('set dlg') # This is used later to get and set

# Get the current command in the gui we just made:
currentcommand = tk.eval ("dialog_get_command " + dlg)

# Set the command in the dialog we just made:
tk.eval ("dialog_set_command " + dlg + " {d.vect map=roads}")
</source>

== FAQ ==

* '''Q:''' Error message "execl() failed: Permission denied" - what to do?
: '''A:''' Be sure that the execute bit of the script is set.

== Links ==

=== General guides ===

* [http://en.wikibooks.org/wiki/Python_Programming/ Wikibook Python Programming]
* [http://www.poromenos.org/tutorials/python Quick Python tutorial] for programmers of other languages
*: [http://wiki.python.org/moin/BeginnersGuide/Programmers More Python tutorials] for programmers
* [http://www.python.org/dev/peps/pep-0008/ Python programming style guide]
* [http://wiki.python.org/moin/PythonEditors Python Editors]

=== Programming ===

* Python and GRASS:
** Library interfaces: [GRASS Python Scripting Library http://grass.osgeo.org/programming7/pythonlib.html]
** Graphical user interface (GIU): [GRASS wxPython-based GUI http://grass.osgeo.org/programming7/wxpythonlib.html]
** PyWPS, GRASS-Web Processing Service: [[WPS]]

* Python and OSGeo:
** [http://wiki.osgeo.org/wiki/OSGeo_Python_Library OSGeo Python Library]

* Python and GDAL/OGR:
** [http://mapserver.gis.umn.edu/community/conferences/MUM3/workshop/python Open Source Python GIS Hacks Mum'03]
** http://hobu.biz/software/OSGIS_Hacks - Python OSGIS Hacks '05
** http://zcologia.com/news/categorylist_html?cat_id=8
** http://www.perrygeo.net/wordpress/?p=4

* Python bindings to PROJ:
** http://www.cdc.noaa.gov/people/jeffrey.s.whitaker/python/pyproj.html

* Python and GIS:
** [http://gispython.org/ Open Source GIS-Python Laboratory]

* Python and Statistics:
** [http://rpy.sourceforge.net/ RPy] - Python interface to the R-statistics programming language

* Bindings:
** SIP (C/C++ bindings generator) http://directory.fsf.org/all/Python-SIP.html
** [http://www.cython.org/ Cython] - C-Extensions for Python (compile where speed is needed)

* Other external projects
** [http://www.scipy.org Scientific Python]
** [http://wiki.python.org/moin/NumericAndScientific Numeric and Scientific]
** [http://w3.pppl.gov/~hammett/comp/python/python.html Info on Python for Scientific Applications]

=== Presentations ===

* http://www.slideshare.net/search/slideshow?q=pygrass

From FOSS4G2006:
* [http://www.foss4g2006.org/materialDisplay.py?contribId=136&sessionId=48&materialId=slides&confId=1 A Python sweeps in the GRASS] - A. Frigeri 2006
* [http://www.foss4g2006.org/materialDisplay.py?contribId=67&sessionId=48&materialId=slides&confId=1 GRASS goes web: PyWPS] - J. Cepicky 2006 (see also [[WPS]])

=== References ===

* Zambelli, P., Gebbert, S., Ciolli, M., 2013. ''Pygrass: An Object Oriented Python Application Programming Interface (API) for Geographic Resources Analysis Support System (GRASS) Geographic Information System (GIS)''. ISPRS International Journal of Geo-Information 2, 201–219. ([http://dx.doi.org/10.3390/ijgi2010201 DOI] | [http://www.mdpi.com/2220-9964/2/1/201/pdf PDF])

=== Python SIGs ===
Python Special Interest Groups are focused collaborative efforts to develop, improve, or maintain specific Python resources. Each SIG has a charter, a coordinator, a mailing list, and a directory on the Python website. SIG membership is informal, defined by subscription to the SIG's mailing list. Anyone can join a SIG, and participate in the development discussions via the SIG's mailing list. Below is the list of currently active Python SIGs, with links to their resources.

See more at http://www.python.org/community/sigs/

{{Python}}

Atmospheric correction

2016-01-19T11:37:39Z

⚠️Andbal: added grass7.x code

See also the [[LANDSAT]] wiki page.

'''Aim of this tutorial:'''

* Transforming a Landsat ETM+ image from DN values to radiance values
* Conducting an atmospheric correction of a Landsat image using the {{cmd|i.atcorr}} module in GRASS GIS and convert radiance to reflectance values

== Introduction ==

Standard steps for pre-processing satellite images are:

# Converting '''Digital Numbers''' (DN) to '''Top-of-Atmosphere Radiances''' (ToAR)
#* for Landsat imagery, apply the equation below or use the {{cmd|i.landsat.toar}} module
#* DNs -- which are actually the pixel values -- are the result of the quantified amount of energy observed and measured at the sensor
#* the {{cmd|i.landsat.toar}} module supports all Landsat sensors icluding MSS, TM and ETM+
# Converting '''ToA Radiances''' to '''Top-of-Canopy Reflectances''' (ToCR) → use the {{cmd|i.atcorr}} module
#* the conversion to ''Reflectances'' attempts to assess for and remove atmospheric effects, a process known as ''atmospheric correction''
#* the ''atmospheric correction'' step is quite a common topic of discussion in Remote Sensing. However, depending on the quality of the acquisition(s), it can be a complex process that does not always offer a great benefit, in terms of accuracy for subsequent processing steps (e.g. image segmentation and classification).

'''Remember''' to check if water areas are represented by values <code>>0</code>, since reflectance is always positive. If you encounter negative reflectance values, you have a "cornichon" (pickle in English ;-0). This means that the transformation equations used, do not correspond to the image processed.

An overview of options to get "corrected" Spectral Reflectances:

<nowiki>+------------------------------------------------------------------------------+
| |
| Digital Numbers |
| | |
| +-----v-----+ |
| | i.*.toar | ---> Reflectance |
| +-----+-----+ (uncorrected)~~~(DOS methods)--+ |
| | + | |
| (-r flag) (-r flag) +--> "Corrected" |
| | | | Reflectance |
| v +-----v----+ | |
| Radiance ------> i.atcorr +-------------------+ |
| +----------+ |
| |
+------------------------------------------------------------------------------+</nowiki>

== Requirements ==

* GRASS 6.4.0 or higher
* North Carolina sample dataset (location): https://grass.osgeo.org/download/sample-data/

After downloading of the North Carolina dataset, it must be unzipped into your GIS database. When starting GRASS GIS, select this GIS database, open the ‘nc_spm_08’ dataset as LOCATION, and ‘PERMANENT’ as MAPSET.

== Calculating radiance values ==

The NC data set mapset contains, amongst others, a Landsat ETM+ image of 24th of May 2002. Every pixel of this image contains a DN or grey value. In order to be able to make calculations with satellite imagery, or compare values amongst different sensors, these values have to be converted to radiances or reflectances. The formulas to do this conversion are described here for Landsat images (or use {{cmd|i.landsat.toar}}).

The conversion from DN values to top-of-atmosphere radiances is done as follows:

<math>L\lambda = \frac{L_{MAX}λ - L_{MIN}λ}{QCAL_{MAX} - QCAL_{MIN}} * (QCAL - QCAL_{MIN}) + L_{MIN}\lambda</math>

<ul>where:

* <math>Lλ</math> - spectral Radiance at the sensor's aperture in watts/(meter squared * ster * μm), the apparent radiance as seen by the satellite sensor ;
* <math>QCAL</math> - the quantized calibrated pixel value in <math>DN</math> ;
* <math>LMINλ</math> - the spectral radiance that is scaled to <math>QCALMIN</math> in watts/(meter squared * ster * μm) ;
* <math>LMAXλ</math> - the spectral radiance that is scaled to <math>QCALMAX</math> in watts/(meter squared * ster * μm) ;
* <math>QCALMINM</math> - the minimum quantized calibrated pixel value (corresponding to <math>LMINλ</math>) in <math>DN</math> ;
* <math>QCALMAX</math> - the maximum quantized calibrated pixel value (corresponding to <math>LMAXλ</math>) in <math>DN=255</math>.
</ul>

<math>LMINλ</math> and <math>LMAXλ</math> are the radiances related to the minimal and maximal <math>DN</math> value, and are reported in the metadata file for each image, or in the table 1. High gain or low gain is also reported in the metadata file of each Landsat image. The minimal <math>DN</math> value (<math>QCALMIN</math>) is 1 for Landsat ETM+ images (1), and the maximal <math>DN</math> value (<math>QCALMAX</math>) is 255. <math>QCAL</math> is the <math>DN</math> value for every separate pixel in the Landsat image.

Accessing the metadata:

{{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_xx</source>

Under ‘comments’, the maximal and minimal radiance (LMAX and LMIN) for each band are given.

Power users may use this for nicely readable output (example for channel 1):

<source lang="bash" enclose=none>CHAN=1</source>
{{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_${CHAN}0 -h | tr '\n' ' ' | sed 's+ ++g' | tr ':' '\n' | grep "LMIN_BAND${CHAN}\|LMAX_BAND${CHAN}"
</source>

Conversion to radiance (this calculation is done for band 1, for the other bands, the numbers in italics need to be replaced with their correct values):

{{cmd|g.region}} <source lang="bash" enclose=none>rast=lsat7_2002_10 -p</source>
{{cmd|r.mapcalc}} <source lang="bash" enclose=none>"lsat7_2002_10_rad=((''191.6'' - (''-6.2'')) / (255.0 - 1.0)) * (lsat7_2002_10 - 1.0) + (''-6.2'')"</source>

For faster computations, use an integer DEM:

{{cmd|r.mapcalc}} <source lang="bash" enclose=none>"elev_int = round(elevation)"</source>

Find mean elevation to initialize computations (needed in 'icnd' file below)

{{cmd|r.univar}} <source lang="bash" enclose=none>elev_int</source>

== Estimating the overpass time ==
This value is needed for the control file. For details, refer to Wikipedia's entry on [http://en.wikipedia.org/wiki/Decimal_time#Scientific_decimal_time Scientific decimal time].
In case this parameter is not reported in the metadata file, we have two ways to obtain it.

=== From sun position ===

The satellite overpass time can be estimated rather precisely from the sun position reported in metadata using {{cmd|r.sunmask}}.
The [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz metadata file] for the following example contains: <source lang="bash" enclose="none">SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999</source>.

# iteratively change hour and minutes to obtain a good result, timezone needs to be correct of course
{{cmd|r.sunmask}} <source lang="bash" enclose=none>-s elev=elevation out=dummy year=2002 month=5 day=24 hour=10 min=42 sec=7 timezone=-5</source>

The above command reports:

Solar position: sun azimuth: 121.342461, sun angle above horz.(refraction corrected): 65.396652

The resulting overpass local time <code>10:42:07</code> corresponds to <code>15:42</code> in GMT which corresponds to <code>15.70</code> in decimal GMT hours (decimal minutes: <math>42 * 100 / 60</math>).

=== From NASA web tool ===

The [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi NASA LaRC Satellite Overpass Predictor] can compute satellite overpass time from d''ate of acquisition'' and ''scene center coordinates''. The [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz metadata file] for the following example contains: <source lang="bash" enclose="none">ACQUISITION_DATE = 2002-05-24, SCENE_CENTER_LAT = +36.0512847, SCENE_CENTER_LON = -79.3280820</source>.

# Select the World zone in the bottom-right side of the page and click <code>Go >></code>
# Insert Lat/Long coordinates (decimal degrees, without ''plus'' sign but with ''minus'' sign)
# Select the proper satellite
# Insert the date of acquisition
# ''Optional:'' select "Day" or "Night" to restrict computation to (local) day/night overpass of the satellite

--- For LANDSAT only ---

USGS provide the [http://landsat.usgs.gov/consumer.php Landsat Bulk Metadata Service] web page where is possible to extract metadata (included overpass time) for all the Landsat missions.

== Atmospheric correction ==

Radiance values can be converted to reflectance values in undergoing an atmospheric correction by applying the <math>6S</math> algorithm. The algorithm will transform the Top-of-Atmosphere radiance values to Top-of-Canopy reflectance values using predefined information on the aerosol content and atmospheric composition of the time the image was taken. What follows describes the method to use this algorithm in GRASS GIS. Again, only calculations for band 1 are demonstrated here, and the number(s) to be changed for indicating other spectral bands, are red. The <code>icnd_lsat1.txt</code> control file consists of the following parameters, and is written with a text editor:

8 # indicates that it is an ETM+ image
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N
2 # the midlatitude summer atmospheric model
1 # the continental aerosol model
50 # the visibility for the aerosol model [km]
-0.110 # the terrain lies 110 meters above sea level [km] * -1
-1000 # image taken of a satellite sensor (1000)
61 # spectral band, here 1: blue

Alternatively, in case visibility maps are not available, an Aerosol Optical Depth estimation at 550nm can be used. The control file should be altered accordingly by entering 0 for the visibility and, in a following line, entering the aerosol optical depth. The control file should be structured as follows:

8 # indicates that it is an ETM+ image
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N
2 # the midlatitude summer atmospheric model
1 # the continental aerosol model
0 # set visibility = 0
0.112 # e.g. aerosol optical depth
-0.110 # the terrain lies 110 meters above sea level [km] * -1
-1000 # image taken of a satellite sensor (1000)
61 # spectral band, here 1: blue

This file is then used in the {{cmd|i.atcorr}} module:

{{cmd|i.atcorr}} <source lang="bash" enclose=none>-a -o iimg=lsat7_2002_10_rad ialt=elev_int icnd=icdn_lsat1.txt oimg=lsat7_2002_10_atcorr</source>

or, using GRASS 7.x:

{{cmd|i.atcorr}} <source lang="bash" enclose=none>-a input=lsat7_2002_10_rad elevation=elev_int parameters=icdn_lsat1.txt output=lsat7_2002_10_atcorr</source>

Where:
* <tt>'''-a'''</tt> refers to a Landsat image taken after July 2000
* <tt>'''-o'''</tt> activates the cache acceleration (GRASS 6.x only)
* <tt>'''iimg/input'''</tt> is the image to be corrected
* <tt>'''ialt/elevation'''</tt> is the altitude map which overrides the initialization value of 110 meters
* <tt>'''icnd/parameters'''</tt> is the path to the icnd.txt file
* <tt>'''oimg/output'''</tt> is the name of the output image

Detailed information can be found in {{cmd|i.atcorr}}'s manual.

'''Note,''' {{cmd|i.atcorr}} wont complain if the defined elevation model (parameter <code>elevation=</code>) does not cover the extent of the Landsat scene that is the subject of the process -- it will result in NaN's. ''[Maybe this can be requested as an enhancement?]''

== Sources for aerosol optical depth (AOD) estimations ==

* [http://aeronet.gsfc.nasa.gov AERONET] - provides globally distributed observations of spectral aerosol optical depth (AOD), inversion products, and precipitable water in diverse aerosol regimes.
* MODIS MOD08 product [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 MOD08 - Gridded Atmospheric Product].
** Daily MOD08 products can be downloaded directly from NASA's FTP server: [ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/ '''ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/'''].
** File and layer specifications are available, for example, at [http://www.atmos.washington.edu/~robwood/modis/filespecs/MOD08_D3.CDL.fs MODIS HDF File Specification MOD08_D3: MODIS Level 3 Daily Atmosphere Gridded Product]
** The layer '''Optical_Depth_Land_And_Ocean_Mean''', described as ''Aerosol Optical Thickness at 0.55 microns for both Ocean (best) and Land (corrected): Mean'', can be probably used as an input for {{cmd|i.atcorr}}. '''Note,''' the layer's values have to be scaled by using the scaling factor ''Optical_Depth_Land_And_Ocean_Mean:scale_factor ='' '''0.001d''' as reported in the specifications.
* [http://www.cesbio.ups-tlse.fr/multitemp/?p=1710 How to estimate Aerosol Optical Thickness]

=== TοDο ===
* understand if the [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=04 MODIS Aerosol Product (MOD04)] could be useful to estimate the visibility for the aerosol model.
* provide an example on using [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 MODIS Gridded Atmospheric Product (MOD08)].
* explain if it is possible and how to use ''Rayleigh Optical Depths estimated at 0km altitude for six different atmosphere models'' from the paper [http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-34-15-2765 Rayleigh-scattering calculations for the terrestrial atmosphere] by Anthony Buchholz, table 4, page 2270

== References ==

=== 6S algorithm ===

* [http://6s.ltdri.org 6S Web site]
* [http://modis-sr.ltdri.org/ Land Surface Reflectance Science Computing Facility website - 6S]

=== Spectral radiance/reflectance ===

* [http://landsat.usgs.gov/how_is_radiance_calculated.php How is radiance calculated?], question on USGS' [http://landsat.usgs.gov/tools_faq.php Landsat Missions - Frequently Asked Questions]
* (1) http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html
* (2) http://landsat.usgs.gov/documents/L5TM_postcal_v11.pdf

=== Whitepapers ===

* Various [http://apollomapping.com/about-us/whitepapers Whitepapers] about conversion and atmospheric correction of High Resolution Satellite Imagery
* Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P., Macomber, S.A., 2001. Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects? Remote Sensing of Environment 75, 230-244. [http://www.unc.edu/courses/2008spring/geog/577/001/www/Song01_RSE.pdf PDF]

=== Related on-line tools ===

* [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi NASA LaRC Satellite Overpass Predictor]
* [http://landsat.usgs.gov/consumer.php Landsat Bulk Metadata Service]
* NASA's [http://atmcorr.gsfc.nasa.gov/ Atmospheric Correction Parameter Calculator]

=== How to add new sensors to i.atcorr ===

How to add new sensors to {{cmd|i.atcorr}}:
<ul>
* see http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/README
</ul>

[[Category: Documentation]]
[[Category: Tutorial]]
[[Category: Image processing]]

Atmospheric correction

2015-12-11T10:43:17Z

⚠️Andbal: Corrected link and update of overpass time sub-section

See also the [[LANDSAT]] wiki page.

'''Aim of this tutorial:'''

* Transforming a Landsat ETM+ image from DN values to radiance values
* Conducting an atmospheric correction of a Landsat image using the {{cmd|i.atcorr}} module in GRASS GIS and convert radiance to reflectance values

== Introduction ==

Standard steps for pre-processing satellite images are:

# Converting '''Digital Numbers''' (DN) to '''Top-of-Atmosphere Radiances''' (ToAR)
#* for Landsat imagery, apply the equation below or use the {{cmd|i.landsat.toar}} module
#* DNs -- which are actually the pixel values -- are the result of the quantified amount of energy observed and measured at the sensor
#* the {{cmd|i.landsat.toar}} module supports all Landsat sensors icluding MSS, TM and ETM+
# Converting '''ToA Radiances''' to '''Top-of-Canopy Reflectances''' (ToCR) → use the {{cmd|i.atcorr}} module
#* the conversion to ''Reflectances'' attempts to assess for and remove atmospheric effects, a process known as ''atmospheric correction''
#* the ''atmospheric correction'' step is quite a common topic of discussion in Remote Sensing. However, depending on the quality of the acquisition(s), it can be a complex process that does not always offer a great benefit, in terms of accuracy for subsequent processing steps (e.g. image segmentation and classification).

'''Remember''' to check if water areas are represented by values <code>>0</code>, since reflectance is always positive. If you encounter negative reflectance values, you have a "cornichon" (pickle in English ;-0). This means that the transformation equations used, do not correspond to the image processed.

An overview of options to get "corrected" Spectral Reflectances:

<nowiki>+------------------------------------------------------------------------------+
| |
| Digital Numbers |
| | |
| +-----v-----+ |
| | i.*.toar | ---> Reflectance |
| +-----+-----+ (uncorrected)~~~(DOS methods)--+ |
| | + | |
| (-r flag) (-r flag) +--> "Corrected" |
| | | | Reflectance |
| v +-----v----+ | |
| Radiance ------> i.atcorr +-------------------+ |
| +----------+ |
| |
+------------------------------------------------------------------------------+</nowiki>

== Requirements ==

* GRASS 6.4.0 or higher
* North Carolina sample dataset (location): https://grass.osgeo.org/download/sample-data/

After downloading of the North Carolina dataset, it must be unzipped into your GIS database. When starting GRASS GIS, select this GIS database, open the ‘nc_spm_08’ dataset as LOCATION, and ‘PERMANENT’ as MAPSET.

== Calculating radiance values ==

The NC data set mapset contains, amongst others, a Landsat ETM+ image of 24th of May 2002. Every pixel of this image contains a DN or grey value. In order to be able to make calculations with satellite imagery, or compare values amongst different sensors, these values have to be converted to radiances or reflectances. The formulas to do this conversion are described here for Landsat images (or use {{cmd|i.landsat.toar}}).

The conversion from DN values to top-of-atmosphere radiances is done as follows:

<math>L\lambda = \frac{L_{MAX}λ - L_{MIN}λ}{QCAL_{MAX} - QCAL_{MIN}} * (QCAL - QCAL_{MIN}) + L_{MIN}\lambda</math>

<ul>where:

* <math>Lλ</math> - spectral Radiance at the sensor's aperture in watts/(meter squared * ster * μm), the apparent radiance as seen by the satellite sensor ;
* <math>QCAL</math> - the quantized calibrated pixel value in <math>DN</math> ;
* <math>LMINλ</math> - the spectral radiance that is scaled to <math>QCALMIN</math> in watts/(meter squared * ster * μm) ;
* <math>LMAXλ</math> - the spectral radiance that is scaled to <math>QCALMAX</math> in watts/(meter squared * ster * μm) ;
* <math>QCALMINM</math> - the minimum quantized calibrated pixel value (corresponding to <math>LMINλ</math>) in <math>DN</math> ;
* <math>QCALMAX</math> - the maximum quantized calibrated pixel value (corresponding to <math>LMAXλ</math>) in <math>DN=255</math>.
</ul>

<math>LMINλ</math> and <math>LMAXλ</math> are the radiances related to the minimal and maximal <math>DN</math> value, and are reported in the metadata file for each image, or in the table 1. High gain or low gain is also reported in the metadata file of each Landsat image. The minimal <math>DN</math> value (<math>QCALMIN</math>) is 1 for Landsat ETM+ images (1), and the maximal <math>DN</math> value (<math>QCALMAX</math>) is 255. <math>QCAL</math> is the <math>DN</math> value for every separate pixel in the Landsat image.

Accessing the metadata:

{{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_xx</source>

Under ‘comments’, the maximal and minimal radiance (LMAX and LMIN) for each band are given.

Power users may use this for nicely readable output (example for channel 1):

<source lang="bash" enclose=none>CHAN=1</source>
{{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_${CHAN}0 -h | tr '\n' ' ' | sed 's+ ++g' | tr ':' '\n' | grep "LMIN_BAND${CHAN}\|LMAX_BAND${CHAN}"
</source>

Conversion to radiance (this calculation is done for band 1, for the other bands, the numbers in italics need to be replaced with their correct values):

{{cmd|g.region}} <source lang="bash" enclose=none>rast=lsat7_2002_10 -p</source>
{{cmd|r.mapcalc}} <source lang="bash" enclose=none>"lsat7_2002_10_rad=((''191.6'' - (''-6.2'')) / (255.0 - 1.0)) * (lsat7_2002_10 - 1.0) + (''-6.2'')"</source>

For faster computations, use an integer DEM:

{{cmd|r.mapcalc}} <source lang="bash" enclose=none>"elev_int = round(elevation)"</source>

Find mean elevation to initialize computations (needed in 'icnd' file below)

{{cmd|r.univar}} <source lang="bash" enclose=none>elev_int</source>

== Estimating the overpass time ==
This value is needed for the control file. For details, refer to Wikipedia's entry on [http://en.wikipedia.org/wiki/Decimal_time#Scientific_decimal_time Scientific decimal time].
In case this parameter is not reported in the metadata file, we have two ways to obtain it.

=== From sun position ===

The satellite overpass time can be estimated rather precisely from the sun position reported in metadata using {{cmd|r.sunmask}}.
The [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz metadata file] for the following example contains: <source lang="bash" enclose="none">SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999</source>.

# iteratively change hour and minutes to obtain a good result, timezone needs to be correct of course
{{cmd|r.sunmask}} <source lang="bash" enclose=none>-s elev=elevation out=dummy year=2002 month=5 day=24 hour=10 min=42 sec=7 timezone=-5</source>

The above command reports:

Solar position: sun azimuth: 121.342461, sun angle above horz.(refraction corrected): 65.396652

The resulting overpass local time <code>10:42:07</code> corresponds to <code>15:42</code> in GMT which corresponds to <code>15.70</code> in decimal GMT hours (decimal minutes: <math>42 * 100 / 60</math>).

=== From NASA web tool ===

The [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi NASA LaRC Satellite Overpass Predictor] can compute satellite overpass time from d''ate of acquisition'' and ''scene center coordinates''. The [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz metadata file] for the following example contains: <source lang="bash" enclose="none">ACQUISITION_DATE = 2002-05-24, SCENE_CENTER_LAT = +36.0512847, SCENE_CENTER_LON = -79.3280820</source>.

# Select the World zone in the bottom-right side of the page and click <code>Go >></code>
# Insert Lat/Long coordinates (decimal degrees, without ''plus'' sign but with ''minus'' sign)
# Select the proper satellite
# Insert the date of acquisition
# ''Optional:'' select "Day" or "Night" to restrict computation to (local) day/night overpass of the satellite

== Atmospheric correction ==

Radiance values can be converted to reflectance values in undergoing an atmospheric correction by applying the <math>6S</math> algorithm. The algorithm will transform the Top-of-Atmosphere radiance values to Top-of-Canopy reflectance values using predefined information on the aerosol content and atmospheric composition of the time the image was taken. What follows describes the method to use this algorithm in GRASS GIS. Again, only calculations for band 1 are demonstrated here, and the number(s) to be changed for indicating other spectral bands, are red. The <code>icnd_lsat1.txt</code> control file consists of the following parameters, and is written with a text editor:

8 # indicates that it is an ETM+ image
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N
2 # the midlatitude summer atmospheric model
1 # the continental aerosol model
50 # the visibility for the aerosol model [km]
-0.110 # the terrain lies 110 meters above sea level [km] * -1
-1000 # image taken of a satellite sensor (1000)
61 # spectral band, here 1: blue

Alternatively, in case visibility maps are not available, an Aerosol Optical Depth estimation at 550nm can be used. The control file should be altered accordingly by entering 0 for the visibility and, in a following line, entering the aerosol optical depth. The control file should be structured as follows:

8 # indicates that it is an ETM+ image
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N
2 # the midlatitude summer atmospheric model
1 # the continental aerosol model
0 # set visibility = 0
0.112 # e.g. aerosol optical depth
-0.110 # the terrain lies 110 meters above sea level [km] * -1
-1000 # image taken of a satellite sensor (1000)
61 # spectral band, here 1: blue

This file is then used in the {{cmd|i.atcorr}} module:

{{cmd|i.atcorr}} <source lang="bash" enclose=none>-a -o iimg=lsat7_2002_10_rad ialt=elev_int icnd=icdn_lsat1.txt oimg=lsat7_2002_10_atcorr</source>

Where:
* <tt>'''-a'''</tt> refers to a Landsat image taken after July 2000
* <tt>'''-o'''</tt> activates the cache acceleration
* <tt>'''iimg'''</tt> is the image to be corrected
* <tt>'''ialt'''</tt> is the altitude map which overrides the initialization value of 110 meters
* <tt>'''icnd'''</tt> is the path to the icnd.txt file
* <tt>'''oimg'''</tt> is the name of the output image

Detailed information can be found in {{cmd|i.atcorr}}'s manual.

'''Note,''' {{cmd|i.atcorr}} wont complain if the defined elevation model (parameter <code>elevation=</code>) does not cover the extent of the Landsat scene that is the subject of the process -- it will result in NaN's. ''[Maybe this can be requested as an enhancement?]''

== Sources for aerosol optical depth (AOD) estimations ==

* [http://aeronet.gsfc.nasa.gov AERONET] - provides globally distributed observations of spectral aerosol optical depth (AOD), inversion products, and precipitable water in diverse aerosol regimes.
* MODIS MOD08 product [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 MOD08 - Gridded Atmospheric Product].
** Daily MOD08 products can be downloaded directly from NASA's FTP server: [ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/ '''ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/'''].
** File and layer specifications are available, for example, at [http://www.atmos.washington.edu/~robwood/modis/filespecs/MOD08_D3.CDL.fs MODIS HDF File Specification MOD08_D3: MODIS Level 3 Daily Atmosphere Gridded Product]
** The layer '''Optical_Depth_Land_And_Ocean_Mean''', described as ''Aerosol Optical Thickness at 0.55 microns for both Ocean (best) and Land (corrected): Mean'', can be probably used as an input for {{cmd|i.atcorr}}. '''Note,''' the layer's values have to be scaled by using the scaling factor ''Optical_Depth_Land_And_Ocean_Mean:scale_factor ='' '''0.001d''' as reported in the specifications.
* [http://www.cesbio.ups-tlse.fr/multitemp/?p=1710 How to estimate Aerosol Optical Thickness]

=== TοDο ===
* understand if the [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=04 MODIS Aerosol Product (MOD04)] could be useful to estimate the visibility for the aerosol model.
* provide an example on using [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 MODIS Gridded Atmospheric Product (MOD08)].
* explain if it is possible and how to use ''Rayleigh Optical Depths estimated at 0km altitude for six different atmosphere models'' from the paper [http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-34-15-2765 Rayleigh-scattering calculations for the terrestrial atmosphere] by Anthony Buchholz, table 4, page 2270

== References ==

=== 6S algorithm ===

* [http://6s.ltdri.org 6S Web site]
* [http://modis-sr.ltdri.org/ Land Surface Reflectance Science Computing Facility website - 6S]

=== Spectral radiance/reflectance ===

* [http://landsat.usgs.gov/how_is_radiance_calculated.php How is radiance calculated?], question on USGS' [http://landsat.usgs.gov/tools_faq.php Landsat Missions - Frequently Asked Questions]
* (1) http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html
* (2) http://landsat.usgs.gov/documents/L5TM_postcal_v11.pdf

=== Whitepapers ===

* Various [http://apollomapping.com/about-us/whitepapers Whitepapers] about conversion and atmospheric correction of High Resolution Satellite Imagery
* Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P., Macomber, S.A., 2001. Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects? Remote Sensing of Environment 75, 230-244. [http://www.unc.edu/courses/2008spring/geog/577/001/www/Song01_RSE.pdf PDF]

=== Related on-line tools ===

* [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi NASA LaRC Satellite Overpass Predictor]
* NASA's [http://atmcorr.gsfc.nasa.gov/ Atmospheric Correction Parameter Calculator]

=== How to add new sensors to i.atcorr ===

How to add new sensors to {{cmd|i.atcorr}}:
<ul>
* see http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/README
</ul>

[[Category: Documentation]]
[[Category: Tutorial]]
[[Category: Image processing]]

LANDSAT

2015-12-11T09:21:13Z

⚠️Andbal: Corrected link to manual page and external site

== Data Availability==

Landsat imagery can be obtained from [http://landsat.usgs.gov/ USGS' LANDSAT archive]:

* [http://glovis.usgs.gov USGS GloVis]
* [http://earthexplorer.usgs.gov/ EarthExplorer]
* [http://earthexplorer.usgs.gov/bulk/help Bulk Download Orders] using the [http://earthexplorer.usgs.gov/bulk/ Bulk Download Application]
* [http://landsatlook.usgs.gov/ LandsatLook Viewer]

See also the [[Global datasets]] wiki page

On 06/12/13, USGS [http://landsat.usgs.gov/about_LU_Vol_7_Issue_4.php announced] the availability of [http://landsat.usgs.gov/landsat8.php LANDSAT 8]. See also [http://landsat.usgs.gov/band_designations_landsat_satellites.php this page].

== Modules overview ==

GRASS GIS features a complete set of native modules and various add-ons for pre- and post-processing Landsat satellite imagery. The following lists offer an overview of related modules and add-ons.

=== Generic modules applicable to Landsat ===

* {{cmd|d.rgb}} → display 3-band data
* {{cmd|r.composite}} → flatten 3-bands of data into a single image (lossy, maybe used in combination with {{cmd|i.landsat.rgb|version=64}} (GRASS 6.x) | {{cmd|i.colors.enhance|version=70}} (GRASS 7.x)
* {{cmd|i.atcorr}} → correct top of atmosphere to surface reflectance - see also the [[Atmospheric correction]] wiki page
* {{cmd|i.topo.corr}} → used to topographically correct reflectance from imagery files, e.g. obtained with {{cmd|i.landsat.toar}}, using a sun illumination terrain model
* {{cmd|i.oif}} → calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
* {{AddonCmd|i.histo.match}} (addon) → histogram matching of two or more raster maps (in grass 7) ''Note, the module works with integer values and does not accept the "." character as part of the raster map's name!''

=== Landsat specific modules ===

* {{cmd|i.landsat.rgb|version=64}} (GRASS 6.x) | {{cmd|i.colors.enhance|version=70}} (GRASS 7.x) → auto-enhance colors
* {{cmd|i.landsat.toar}} → convert DN to top of atmosphere radiance
* {{cmd|i.landsat.acca}} → cloud cover assessment

=== Landsat specific GRASS AddOns ===

* {{AddonCmd|i.landsat.trim|version=64}} → trims border fringes for each band separately or with the MASK where coverage exists for all bands
* {{AddonCmd|i.landsat.dehaze|version=64}} (addon) → haze removal

== Pre-Processing ==

=== Overview ===

Typically, pre-processing Landsat imagery comprises the following steps:

# '''import''' in the database → {{cmd|r.in.gdal}}
# geometrically & orthometrically correct imagery
#* already done for L1T products, read more at [http://landsat.usgs.gov//Landsat_Processing_Details.php USGS' Landsat Information Products] webpage
# optionally, automatically '''cut-off border fringes''' → {{AddonCmd|i.landsat.trim|version=64}}
#* of course one can use the official WRS2 Path/Row vector tiles to manually trim border fringes
# optionally, '''denoise''' for obvious/intensive salt & pepper effects, stripes, etc.
#* for example by applying Principal Components Analysis as a denoising technique → {{cmd|i.pca}} ''<<< Re-order this step?''
# '''convert''' the '''Digital Numbers''' (DNs) to '''Top-of-Atmosphere Radiances/Reflectances''' (ToARs) → {{cmd|i.landsat.toar}}
# optionally, '''correct for atmospheric effects''' → {{cmd|i.atcorr}}
#* that is, accounting for distorting atmospheric effects and estimating actual reflectances as they would have been measured on the ground
#* also described as conversion to Top-of-Canopy Reflectances (ToCRs)
# '''assessing cloud cover''' → {{cmd|i.landsat.acca}}
#* optionally, detect and remove clouds shadows as well
# '''topographically normalise''' imagery → {{cmd|i.topo.corr}}
#* also known as topographic correction, that is, accounting for illumination differences due to the acquisition's geometry
# '''radiometrically normalise''' → one approach via {{AddonCmd|i.histo.match}} (in '''grass 7''', addon), also known as relative radiometric normalisation -- one approach is the ''histogram matching'' technique of two or more raster maps
#* '''ToDo/More techniques?'''

See also
* [http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ Processing Landsat 8 data in GRASS GIS 7: Import and visualization]

=== Importing data ===

Importing Landsat spectral bands in GRASS GIS' data base can be done both from the Graphical User Interface or from the command line

# Open GRASS GIS, select ''Location Wizard'' in order to create a new location from georeferenced file
# Use the ''Import file tool'', or the {{cmd|r.in.gdal}} module, to import the GeoTIFF files into GRASS GIS.

See also [http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import Automated data import] below.

==== Notes ====

# Most Landsat scenes are delivered in ''north-is-up'' orientation, hence the import process is straightforward.
# If you get <source lang="bash" enclose="none">ERROR: Input map is rotated - cannot import.</source>, the image must be first rotated externally, applying the rotation info stored in the metadata field of the raster image file. For example, the <code>gdalwarp</code> software can be used to transform the map to North-up (note, there are several <code>gdalwarp</code> parameters to select the resampling algorithm):

<ul>
<source lang="bash">
gdalwarp rotated.tif northup.tif
</source>
</ul>

==== Hint: Minimal disk space copies ====

Here's a little trick using {{cmd|r.reclass}} to rename maps (for example, from ''L71074092_09220040924_B'''10''''' to ''L71074092_09220040924_B'''1''''') without touching the data or wasting disk space:

<source lang="bash">
BASE=L71074092_09220040924

for BAND in 10 20 30 40 50 61 70 80; do
BAND1st=`echo $BAND | sed -e 's/0$//'`
r.reclass in="${BASE}_B$BAND" out=$BASE.$BAND1st << EOF
* = *
EOF
done
</source>

==== Automated data import ====

The following ''Python'' script imports Landsat imagery into GRASS' data base. Specifically, the script

* creates an independent Mapset for each Landsat scene
* imports and renames bands of a scene as <tt>B<id></tt>, e.g. B10, B20, ..., B80.
* additionaly sets up the timestamp based on MTL metadata file

'''Note,''' the (newest) official naming pattern for Landsat scenes -- explained in [https://lta.cr.usgs.gov/landsat_dictionary.html USGS' Landsat Data Dictionary] as the [https://lta.cr.usgs.gov/landsat_dictionary.html#entity_id Landsat Scene Identifier] -- and all individual bands that compose a scene -- bands have a suffix which is like _B10, _B20, _B30, etc. -- differ from what some Landsat specific modules expect. For example, the modules {{cmd|i.landsat.toar}} and {{cmd|i.landsat.acca}} expect the bands to follow a naming pattern such as "scenename.1, .2, .3", etc.

To use the script save it as <source lang="bash" enclose="none">import_landsat.py</source> file and make sure it is granted the execution permission.

<source lang=python>
#!/usr/bin/python

import os
import sys
import glob
import grass.script as grass

def get_timestamp(mapset):
try:
metafile = glob.glob(mapset + '/*MTL.txt')[0]
except IndexError:
return

result = dict()
try:
fd = open(metafile)
for line in fd.readlines():
line = line.rstrip('\n')
if len(line) == 0:
continue
if any(x in line for x in ('DATE_ACQUIRED', 'ACQUISITION_DATE')):
result['date'] = line.strip().split('=')[1].strip()
finally:
fd.close()

return result

def import_tifs(mapset):
meta = get_timestamp(mapset)
for file in os.listdir(mapset):
if os.path.splitext(file)[-1] != '.TIF':
continue
ffile = os.path.join(mapset, file)
if ('VCID') in ffile:
name = "".join((os.path.splitext(file)[0].split('_'))[1::2])
else:
name = os.path.splitext(file)[0].split('_')[-1]
if len(name) == 3 and name[-1] == '0':
band = int(name[1:2])
elif len(name) == 3 and name[-1] != '0':
band = int(name[1:3])
else:
band = int(name[-1:])
grass.message('Importing %s -> %s@%s...' % (file, name, mapset))
grass.run_command('g.mapset',
flags = 'c',
mapset = mapset,
quiet = True,
stderr = open(os.devnull, 'w'))
grass.run_command('r.in.gdal',
input = ffile,
output = name,
quiet = True,
overwrite = True,
title = 'band %d' % band)
if meta:
year, month, day = meta['date'].split('-')
if month == '01':
month = 'jan'
elif month == '02':
month = 'feb'
elif month == '03':
month = 'mar'
elif month == '04':
month = 'apr'
elif month == '05':
month = 'may'
elif month == '06':
month = 'jun'
elif month == '07':
month = 'jul'
elif month == '08':
month = 'aug'
elif month == '09':
month = 'sep'
elif month == '10':
month = 'oct'
elif month == '11':
month = 'nov'
elif month == '12':
month = 'dec'

grass.run_command('r.timestamp',
map = name,
date = ' '.join((day, month, year)))

def main():
if len(sys.argv) == 1:
for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
import_tifs(directory)
else:
import_tifs(sys.argv[1])

if __name__ == "__main__":
main()
</source>

'''Example of usage'''

After having collected the Landsat scenes of interest in one directory, the script can be used as follows:

* <source lang="python" enclose="none">./import_landsat.py</source> → the script will walk through a ''pool'' directory that contains unique Landsat scene directories (e.g. three directories named after the official naming pattern: LT51800342011158MOR00 LT51810352009079MTI00 LT51820352009326MTI00) and import all bands of each individual scene in their own Mapset
* <source lang="python" enclose="none">./import_landsat.py LM41890261983200FFF03</source> → the scrip will import only bands of the specified Landsat scene directory

== Post-Processing ==

=== Natural color composites ===

Creating natural (also known as true-) color composites, can be done by

a) using {{cmd|i.landsat.rgb|version=64}} (GRASS 6.x) | {{cmd|i.colors.enhance|version=70}} (GRASS 7.x) to automatically balance the colors for the Red, Green and Blue bands

or

b) equalizing the colors on each of the Red, Green and Blue bands

<ul>
<source lang="bash">
r.colors -e map=band1 color=grey
</source>
</ul>

'''Composites''' can then be produced with the {{cmd|r.composite}} module.

'''Resetting the color tables''' of all bands back to normal greyscale can be done with a for loop:

<source lang="bash">
BASE=L71074092_09220040924

for map in `g.mlist pat="$BASE.[0-8]*"` ; do
r.colors $map color=grey255
done
</source>

See also:
* [http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening]

=== Create a MASK to only show data where coverage exists for all bands ===

<source lang="bash">
BASE=L71074092_09220040924

g.region rast=$BASE.1
r.series in=`g.mlist pat="$BASE.[0-8]*" sep=,` -n out=$BASE.thresh method=threshold thresh=1
r.mapcalc "$BASE.mask = if(isnull($BASE.thresh))"
g.remove $BASE.thresh
</source>

=== Calculate Top-of-Atmosphere Reflectance and band-6 Temperature ===

Calculate ''Top-of-Atmosphere'' reflectance and band-6 temperature using the {{cmd|i.landsat.toar}} module. For details, refer to USGS' [http://landsat.usgs.gov/Landsat_Metadata_Changes.php Landsat Filename and Metadata Changes] dedicated page. For Landsat 8, see also [http://landsat.usgs.gov/about_LU_Vol_7_Issue_4.php#3a here].
<source lang="bash">
i.landsat.toar input_prefix=$BASE output_prefix=${BASE}_toar metfile=${BASE}_MTL.txt
</source>

Note, the resulting temperature map is in Kelvin:
<source lang="bash">
# convert to degree Celsius
r.mapcalc "$BASE.temp_celsius = ${BASE}_toar.6 - 273.15"
r.info -r $BASE.temp_celsius
</source>

=== Haze removal ===

Simple haze removal can be done with {{AddonCmd|i.landsat.dehaze|version=64}}. i.landsat.dehaze applies a bandwise haze correction using tasscap4 (haze) and linear regression.

=== Atmospheric correction ===

See [[Atmospheric correction]]

=== Cloud identification ===

Identify clouds in the image with {{cmd|i.landsat.acca}}:
<source lang="bash">i.landsat.acca -f input_prefix=226_62_toar. output=226_62.acca_cloudmask</source>

Mask out the clouds:
<source lang="bash">r.mapcalc "MASK = if(isnull($BASE.acca_cloudmask))"</source>

== Download sample data ==

=== Preprocessed Landsat-7 data for North Carolina ===

The [[Sample_datasets|North Carolina 2008 sample dataset]] comes with 3 different Landsat scenes:
: (''Wake County -- path: 16 row: 35 for various dates'')

The above import efforts are not needed since the data are already in a GRASS location.

=== Landsat-5: Oct 14, 1987 ===

;Glovis download: LT50160351987287XXX08
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LT50160351987287XXX08&dataset_name=LANDSAT_TM



Provided metadata:
IMAGE_ID=P016R35_5T871014
PATH=16
ROW=35
DATE=10/14/87
PLATFORM=LANDSAT5

=== Landsat-7: Mar 31, 2000 ===

;Glovis download: LE70160352000091EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352000091EDC00&dataset_name=LANDSAT_ETM

Values in the metadata below indicate that the version provided with the NC 2008 dataset's production date was after July 1, 2000, and that the channel gains were <tt>HHHLHLHHL</tt>.

Convert DNs to radiance/temperatures:<source lang="bash">
GRASS> i.landsat.toar -v band=lsat7_2000 sensor=7 date=2000-03-31 \
product_date=2000-07-02 solar_elevation=51.5246529 gain=HHHLHLHHL</source>

Provided metadata:
SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
'''ACQUISITION_DATE=2000-03-31'''
WRS_PATH=16
CPF_FILE_NAME=L7CPF20000101_20000331_12
SUN_AZIMUTH=139.6033279
'''SUN_ELEVATION=51.5246529'''
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

=== Landsat-7: May 24, 2002 ===

;Glovis download: LE70160352002144EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352002144EDC00&dataset_name=LANDSAT_ETM



Provided metadata: (`p016r035_7x20020524.met`)
SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
'''ACQUISITION_DATE=2002-05-24'''
WRS_PATH=016
WRS_ROW=035
SUN_AZIMUTH=120.8810347
'''SUN_ELEVATION=64.7730999'''
QA_PERCENT_MISSING_DATA=0
CLOUD_COVER=0
CPF_FILE_NAME=L7CPF20020401_20020630_03
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

== LANDSAT Image classification ==

See [[Image classification]]

== Time series analysis ==

See [[Time series]]

== Sources ==

* [http://landsat.usgs.gov/tools_access_all_faqs.php Frequently Asked Questions about the Landsat Missions]
* [https://landsat8portal.eo.esa.int/portal/ European LANDSAT 8 data in near-realtime]
* [http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ Processing Landsat 8 data in GRASS GIS 7: Import and visualization]
* [http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening]

[[Category: Documentation]]
[[Category: FAQ]]
[[Category: Landsat]]
[[Category: Image processing]]
[[Category: Import]]
[[Category: Digital Numbers]]
[[Category: Radiance]]
[[Category: Reflectance]]
[[Category: Geodata]]

Convert degree

2015-12-01T10:39:25Z

⚠️Andbal: restored link to external site

== Converting between sexagesimal and decimal degree notations ==

=== Converting from sexagesimal (DMS) to decimal degree notation ===

To convert from sexagesimal degree to decimal degree notation the degree
value is kept, the minutes value is divided by 60, the seconds value is divided
by 3600 and the results are summed. An example of a conversion of 12 degree
45 min 50 sec north (usually denotes as 12:45:50N):

12 degree -> 12.00000 degree
45 min / 60 -> 0.75000 degree
50 sec / 3600 -> 0.01389 degree
---------------- addition of values ----------
12.76389 degree (decimal degree notation)

=== Converting from decimal to sexagesimal degree notation ===

To convert from decimal degree to sexagesimal degree notation the calculation is
done as illustrated by an example for a conversion of 12.7689 degree (decimal
notation):

12.76389 -> 12 degree, remainder 0.76389
0.76389 * 60 -> 45 minutes, remainder 0.8334
0.8334 * 60 -> 50 seconds
----------------------------------------------
12:45:50N degree (sexagesimal degree notation)

The N for north has to be added to complete the sexagesimal notation, since
the value is positive in decimal notation. GRASS supports both notations for
latitude-longitude coordinate system.

== See also ==
* https://www.fcc.gov/encyclopedia/degrees-minutes-seconds-tofrom-decimal-degrees
* [[AddOns#Miscellaneous_Add-ons|utm_which_zone.sh]], a shell script to determine UTM zone from Lat/Lon input

[[Category: FAQ]]