IKONOS: Difference between revisions
m (→Automatising Conversions: script constants, band parameters on script's top) |
m (→Importing data: looping r.in.gdal over GeoTIFF bands) |
||
Line 50: | Line 50: | ||
''ToAdd'' | ''ToAdd'' | ||
<ul> | |||
* Location creation based on georeferenced data | |||
</ul> | |||
Once inside a Location that is defined by the spatial reference system in which the bands of interest are projected, they can be imported with the {{cdm|r.in.gdal}} module. | |||
For example, GeoTIFF files can be imported by looping {{cdm|r.in.gdal}} over all of them | |||
<source lang="bash"> | |||
for TIF in `echo *.tif`; do r.in.gdal in=${TIF} out=${TIF%%.*}; done | |||
</source> | |||
=== Deriving Physical Quantities === | === Deriving Physical Quantities === |
Revision as of 16:09, 29 July 2013
[This page is under construction]
IKONOS is a commercial earth observation satellite. Details about the sensor are provided at Digital Globe's IKONOS Data Sheet
Availability (Sample Data)
- Search for commercial satellite image providers in the internet.
- The Global Land Cover Facility (GLCF) provides four openly available IKONOS scenes of western Sichuan.
- ISPRS provides a small IKONOS data set, fragments from a Panchromatic image as well as from a Stereo product.
Pre-Processing Overview
Typically, multispectral satellite data are converted into physical quantities such as Radiance or Reflectance before they are subjected in multispectral analysis techniques (image interpretation, band arithmetic, vegetation indices, matrix transformations, etc.). The latter can be differentiated in Top of Atmosphere Reflectance (ToAR) which does not account for atmospheric effects (absorption or scattering) and in Top of Canopy Reflectance (ToCR) which introduces a "correction" for atmospheric effects.
In order to derive Reflectance values, likewise as with remotely sensed data acquired by other sensors, IKONOS raw image digital numbers (DNs) need to be converted to at-sensor spectral Radiance values. At-sensor spectral Radiance values are an important input for the equation to derive Reflectance values. Note, Spectal Radiance is the measure of the quantity of radiation that hits the sensor and typically expressed in , that is watts per unit source area, per unit solid angle, and per unit wavelength.
Converting DNs to at-sensor Radiance can be done by using the following equation:
Converting to Top of Atmosphere Reflectance, also referred to as Planetary Reflectance, can be done by using the following equation:
where:
- - Unitless Planetary Reflectance
- - mathematical constant (3.14159265358)
- spectral Radiance at the sensor's aperture, from equation... ToADD
- - Earth-Sun distance in astronomical units, interpolated values
- - Mean solar exoatmospheric irradiance(s) (W/m2/μm), interpolated values
- - Solar_zenith_angle from the image acquisition's metadata
Modules overview
ToAdd
File Formats & Metadata
ToAdd
Pre-Processing
ToAdd
Importing data
ToAdd
- Location creation based on georeferenced data
Once inside a Location that is defined by the spatial reference system in which the bands of interest are projected, they can be imported with the Template:Cdm module.
For example, GeoTIFF files can be imported by looping Template:Cdm over all of them
for TIF in `echo *.tif`; do r.in.gdal in=${TIF} out=${TIF%%.*}; done
Deriving Physical Quantities
Converting Digital Numbers to Radiance/Reflectance requires knowledge about the sensor's specific spectral band parameters. Those are, as extracted from the document IKONOS Planetary Reflectance and Mean Solar Exoatmospheric Irradiance, by Martin Taylor (see references):
Pan_CalCoef=161 ; Pan_Width=403 ; Pan_Esun=1375.8
Blue_CalCoef=728 ; Blue_Width=71.3 ; Blue_Esun=1930.9
Green_CalCoef=720 ; Green_Width=88.6 ; Green_Esun=1854.8
Red_CalCoef=949 ; Red_Width=65.8 ; Red_Esun=1556.5
NIR_CalCoef=843 ; NIR_Width=95.4 ; NIR_Esun=1156.9
The following examples exemplify the conversion of raw Blue band digital numbers into Radiance and Reflectance.
Spectral Radiance
Converting a Blue Band (Digital Numbers) in to Spectral at-sensor Radiance in the correct units to be further used for the conversion in to unitless Reflectance:
r.mapcalc "Blue_Radiance = ( (10000 * IKONOS_Blue_Band_DN) / (728 * 71.3) )"
Planetary Reflectance
The equation to derive Reflectance values incorporates in addition:
- The mathematical constant Pi
PI=3.14159265358
. - The Earth-Sun distance in astronomical units which depends on the acquisition's day of year (DOY -- also referred to as Julian day, Ordinal_date) and can be retrieved from the following spreadsheet <http://landsathandbook.gsfc.nasa.gov/excel_docs/d.xls>.
- The mean solar exoatmospheric irradiance in . See 3rd column of (interplated) values given above.
- The cosine of the Solar Zenith Angle (SZA) at the time of the acquisition. The SZA can be calculated from its complementary Solar Elevation Angle (SEA) given in the image acquisition's metadata.
In the following example we accept as the acquisition's DOY=166
and SEA=52.78880
. Hence, we get the Earth-Sun distance ESD=1.0157675
and SZA = 37.21120 deg
.
Converting the in-Blue spectral band at-sensor Radiance in to Planerary Reflectance:
PI=3.14159265358; ESD=1.0157675; BAND_Esun=1930.9; SZA=37.21120
r.mapcalc"Blue_Reflectance = ( ${PI} * Blue_Radiance * ${ESD}^2 ) / ( ${BAND_Esun} * cos(${SZA}) )"
Automatising Conversions
The conversion process can be scripted to avoid repeating the same steps for each band separately. In bash, such a script might be as the following example. Note, however, in this example script constants, band parameters and acquisition related metadata are hard-coded!
#!/bin/bash
Pan_CalCoef=161 ; Pan_Width=403 ; Pan_Esun=1375.8
Blue_CalCoef=728 ; Blue_Width=71.3 ; Blue_Esun=1930.9
Green_CalCoef=720 ; Green_Width=88.6 ; Green_Esun=1854.8
Red_CalCoef=949 ; Red_Width=65.8 ; Red_Esun=1556.5
NIR_CalCoef=843 ; NIR_Width=95.4 ; NIR_Esun=1156.9
# set constants, band parameters and acquisition related metadata here!
# Pi, first 11 decimals
PI=3.14159265358
# Acquisition's Day of Year and estimated Earth-Sun Distance
DOY=166; ESD=1.0157675
# Sun Elevation Angle during the acquisition
SEA=52.78880
# bands to process
Spectral_Bands="Pan Blue Green Red NIR"
# loop over all bands
for BAND in ${Spectral_Bands}; do
echo "Processing the ${BAND} spectral band"
# set region
g.region rast=${BAND}_DNs
# set band parameters as variables
eval BAND_CalCoef="${BAND}_CalCoef"
eval BAND_Width="${BAND}_Width";
echo "Band Parameters set to CalCoef=${!BAND_CalCoef}, Bandwidth=${!BAND_Width}"
# conversion to Radiance
r.mapcalc "${BAND}_Radiance = ( ( 10^4 * ${BAND}_DNs ) / ( ${!BAND_CalCoef} * ${!BAND_Width} ) )"
# add info
r.support map=${BAND}_Radiance \
title="" \
units="W / m2 / μm / ster" \
description="At-sensor `echo ${BAND}` band spectral Radiance (W/m2/μm/sr)" \
source1='"IKONOS Planetary Reflectance and Mean Solar Exoatmospheric Irradiance", by Martin Taylor, Geoeye'
# get Sun Zenith Angle
SZA=$(echo "90 - ${SEA}" | bc )
eval BAND_Esun="${BAND}_Esun"; echo "Using Esun=${!BAND_Esun}"
# calculate ToAR
r.mapcalc "${BAND}_ToAR = \
( ${PI} * ${BAND}_Radiance * ${ESD}^2 ) / ( ${!BAND_Esun} * cos(${SZA}) )"
# add some metadata
r.support map=${BAND}_ToAR \
title="echo ${BAND} band (Top of Atmosphere Reflectance)" \
units="Unitless" \
description="Top of Atmosphere `echo ${BAND}` band spectral Reflectance (unitless)" \
source1='"IKONOS Planetary Reflectance and Mean Solar Exoatmospheric Irradiance", by Martin Taylor, Geoeye' \
source2="e.g., the Image Provider!" \
history="PI=3.14159265358; ESD=1.0157675; BAND_Esun=1930.9; SZA=37.21120"
done
ToAdd
Atmospheric correction
Post-Processing
Having beforehand satellite image data expressed in physical quantities (radiance or reflectance) is preferred to follow-up with digital image analysis techniques.
Color composites
ToAdd
Pan Sharpening
Pan-Sharpening / Fusion is the process of merging high-resolution panchromatic and lower resolution multi-spectral imagery. GRASS 7 holds a dedicated pan-sharpening module, i.pansharpen which features three techniques for sharpening, namely the Brovey transformation, the classical IHS method and one that is based on Principal Components Analysis (PCA).
One approach inside GRASS-GIS to get an acceptable color-balanced composite image after Pan-sharpening 11-bit IKONOS spectral bands comprises the following steps
- rescale the 11-bit IKONOS spectral bands to 8-bit ranging in [0, 255] (r.rescale)
- pan-sharpen with any of the featured methods (Brovey, IHS, PCA) (i.pansharpen)
- color-balance by using the i.landsat.rgb module or manually adjusting the color tables of the bands of interest
- create a composite image by using the r.compose module
Example Instructions
i.pansharpen works fine with 8-bit raster maps as an input. If the data to be processed are out of this range, that is out of , they can be rescaled to fit into this range by using GRASS' r.rescale module.
Given an IKONOS set of 11-bit spectral bands (Blue, Green, Red, NIR and Pan) ranging between , and then querying for example the Blue band
r.info Blue_DNs -r
, would return
min=0
max=2047
Rescaling the Blue band to range between `[0, 255]`
r.rescale in=Blue_DNs out=Blue_DNs_255 from=0,2047 to=0,255
The same step applies to both the rest of the multi-spectral bands and the Panchromatic band of interest.
As usual when working with GRASS, it is required to set the region of interest, i.e. g.region rast=Blue_DNs_255
to match the extent of the band(s) or else. The resolution itself is taken care in this particular case by the module and the resulting pan-sharpened raster maps will of the same high(er) resolution as the Panchromatic band.
An example command for an IHS-based Pan-Sharpening action might look like
i.pansharpen pan=Pan_DNs_255 ms1=Blue_DNs_255 ms2=Green_DNs_255 ms3=Red_DNs_255 output=sharptest255 sharpen=ihs
After the process completion, the module outputs
...
The following pan-sharpened output maps have been generated:
sharptest255_red
sharptest255_green
sharptest255_blue
To visualize output, run: g.region -p rast=sharptest255.red
d.rgb r=sharptest255_red g=sharptest255_green b=sharptest255_blue
If desired, combine channels into a single RGB map with 'r.composite'.
Normally it should be enough to re-balance the colors after the pan-sharpening action by using for example the i.landsat.rgb module or manual adjustment of each of the three bands that would compose an RGB image.
i.landsat.rgb r=sharptest255_red g=sharptest255_green b=sharptest255_blue -p
Vegetation Indices
ToAdd
IKONOS Image classification
References / Sources
- http://www.apollomapping.com/wp-content/user_uploads/2011/09/IKONOS_Esun_Calculations.pdf
- Ikonos DN Value Conversion to Planetary Reflectance Values, by David Fleming
- Landsat7 Science Data Users Handbook, Chapter 11, Section 3
- Some short presentation about the DN to Reflectance conversion: Calibrated Landsat Digital Number (DN) to Top of Atmosphere (TOA) Reflectance Conversion, by Richard Irish
See also
- GRASS-Wiki page about Image Processing
ToAdd