QuickBird
Overview
- Wikipedia's QuickBird entry
Scenes are 16.5km x 16.5km.
Data is 11-bit integers stored in a 16-bit integer field. Thus band minimum = 0, band maximum = 2047.
Multi-spectral bands (2.44 - 2.88m pixel res)
band 1 = blue 2 = green 3 = red 4 = Near IR
Panchromatic band (61-72cm pixel res)
From the whitepaper Radiance Conversion of QuickBird Data:
QuickBird products are delivered to the customer as radiometrically corrected image pixels (qPixel,Band).
Radiometric correction includes a dark offset subtraction and a non-uniformity correction (detector-to-detector relative gain). Corrected counts are specific to the QuickBird instrument and therefore QuickBird imagery MUST be converted to spectral radiance before radiometric/spectral analysis or comparison with imagery from other sensors in
a radiometric/spectral manner.
Availability (Sample Data)
- Search for commercial satellite image providers in the internet.
- The Global Land Cover Facility (GLCF) provides some QuickBird acqusitions for education or research purposes only.
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, QuickBird raw image digital numbers (DNs), also referred to as radiometrically corrected image pixels need to be converted to Top of Atmosphere Radiance values. Spectral Radiance values are then converted to Reflectance values based on various acquisition related parameters. 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 QuickBird DNs to Top of Atmosphere Radiance, as described by Keith Krause (2005) in Radiometric Use of QuickBird Imagery, Technical Note, can be done by using the following equation:
- : Top-of-Atmosphere Spectral Radiance image pixels [W-m-2-sr-1-μm-1]
- : the absolute radiometric calibration factor [W-m-2-sr-1-count-1] for a given band
- : radiometrically corrected image pixels [counts or Digital Numbers]
- : the effective bandwidth [μm] for a given band
- Converting to Top of Atmosphere Reflectance, also referred to as Planetary Reflectance, can be done by using the following equation:
- - 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
where:
where:
Modules overview
Outside of GRASS
- gdalinfo
- gdalwarp
Inside GRASS GIS
File Formats & Metadata
- GeoTIFF, NITF 2.1, NITF 2.0
- Image bits / pixel: 8 or 16 bits
Pre-Processing
ToAdd
Importing data
Setting up a Location in GRASS GIS requires to know the Spatial Reference System in which the data of interest are projected. Retrieving the projection parameters of QuickBird imagery, can be done by using the `gdalinfo` program.
ToAdd
- Various ways of Location creation based on georeferenced data
Once inside an appropriate Location, importing the data is done with the r.in.gdal module.
For example, GeoTIFF files/bands, can be imported in one go by looping r.in.gdal over all of them
for TIF in `echo *.tif`; do r.in.gdal in=${TIF} out=${TIF%%.*}; done
Cleanup
Rename the NIR band
g.rename QBird_Multichrom.4,QBird_Multichrom.NIR
Red and blue bands are swapped, so fix that by renaming them. (see [1])
g.rename QBird_Multichrom.blue,QBird_Multichrom.realred g.rename QBird_Multichrom.red,QBird_Multichrom.blue g.rename QBird_Multichrom.realred,QBird_Multichrom.red
No data is set to 0, so convert that to NULL with r.null:
for BAND in red green blue NIR ; do r.null QBird_Multichrom.$BAND setnull=0 echo $BAND: r.info -r QBird_Multichrom.$BAND echo done
Set color tables appropriately for 11-bit data:
for BAND in red green blue NIR ; do
r.colors QBird_Multichrom.$BAND color=rules << EOF
0 black
2047 white
EOF
done
It is not necessarily needed, but if you wish to rescale to 8-bit data you can use r.rescale as follows:
for BAND in red green blue NIR ; do r.rescale in=QBird_Multichrom.$BAND from=0,2047 \ out=QBird_Multichrom255.$BAND to=0,255 done
Display the 3-band RGB image:
d.rgb r=QBird_Multichrom.red g=QBird_Multichrom.green b=QBird_Multichrom.blue
If the image is excessively dark you might try using the i.landsat.rgb module to auto-balance the colors and redraw:
i.landsat.rgb r=QBird_Multichrom.red g=QBird_Multichrom.green b=QBird_Multichrom.blue d.redraw
Deriving Physical Quantities
Converting Digital Numbers to Radiance/Reflectance requires knowledge about the sensor's specific spectral band parameters. These are, as extracted from the document ...ToAdd:
# Band-Averaged Solar Spectral Irradiance [W/sq.m./μm]
Pan_Esun=1381.79
Blue_Esun=1924.59
Green_Esun=1843.08
Red_Esun=1574.77
NIR_Esun=1113.71
# Revised K Factors & Effective Bandwidths
# 1st column: K Conversion Factors (& Pan TDI Level) for 16-Bit products [W/sq.m./sr/count]
# 2nd column: Effective Bandwidth [μm] per Spectral Band
K16_Pan10=0.08381880 ; Pan10_Width=0.398
K16_Pan13=0.06447600 ; Pan13_Width=0.398
K16_Pan18=0.04656600 ; Pan18_Width=0.398
K16_Pan24=0.03494440 ; Pan24_Width=0.398
K16_Pan32=0.02618840 ; Pan32_Width=0.398
K16_Blue=0.01604120 ; Blue_Width=0.068
K16_Green=0.01438470 ; Green_Width=0.099
K16_Red=0.01267350 ; Red_Width=0.071
K16_NIR=0.01542420 ; NIR_Width=0.114
# 1st column: k′ Conversion Factors for 8-Bit products
# 2nd column being identical (as above) not repeated!
K8_Pan10=1.02681367
K8_Pan13=1.02848939
K8_Pan18=1.02794702
K8_Pan24=1.02989685
K8_Pan32=1.02739898
K8_Blue=1.12097834
K8_Green=1.37652632
K8_Red=1.30924587
K8_NIR=0.98368622
Attention for 32-bit floating point calculations is requried. As stated, in the same document,
NOTE: conversion equations are to be performed on all pixels in a given band of a QuickBird image and should use 32-bit floating point calculations.
Spectral Radiance
Converting a Blue Band (Digital Numbers) in to Spectral Radiance in the correct units to be further used for the conversion in to unitless Reflectance:
# set the region
g.region rast=Blue_DNs -p
# convert DNs to spectral Radiance values
r.mapcalc "Blue_Radiance = '''ToAdd''' "
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=...
and SEA=...
. Hence, we get the Earth-Sun distance ESD=...
and SZA = ... deg
.
Converting the in-Blue spectral band at-sensor Radiance in to Planerary Reflectance:
PI=3.14159265358; ESD=...; BAND_Esun=...; SZA=...
r.mapcalc"Blue_Reflectance = ( ${PI} * Blue_Radiance * ${ESD}^2 ) / ( ${BAND_Esun} * cos(${SZA}) )"
Automatising Conversions
Scripting the conversion process avoids repeating the same steps for each band separately. Note, however, in the following example, constants, band parameters and acquisition related metadata are hard-coded! Reviewing the code and altering appropriately is required, i.e. checking for the parameters ESD
, SEA
, PanTDI
, K_BAND
.
#!bin/bash
# Pi, first 11 decimals
PI=3.14159265358
# HardCoded MetaData!
# Acquisition's Day of Year and estimated Earth-Sun Distance
DOY=274; ESD=1.001190
# Sun Zenith Angle based on the acquisition's Sun Elevation Angle
SEA=67.8; SZA=$(echo "90 - ${SEA}" | bc )
# which Pan TDI level?
PanTDI="10" # evaluate below...
# Spectral Irradiance [W-m-2-μm-1]
Pan_Esun=1381.79
Blue_Esun=1924.59
Green_Esun=1843.08
Red_Esun=1574.77
NIR_Esun=1113.71
# 1st column: K Conversion Factors for 16-Bit products
# 2nd column: Effective Bandwidth [μm] per Spectral Band
K16_Pan10=0.08381880; Pan10_Width=0.398
K16_Pan13=0.06447600; Pan13_Width=0.398
K16_Pan18=0.04656600; Pan18_Width=0.398
K16_Pan24=0.03494440; Pan24_Width=0.398
K16_Pan32=0.02618840; Pan32_Width=0.398
K16_Blue=0.01604120 ; Blue_Width=0.068
K16_Green=0.01438470; Green_Width=0.099
K16_Red=0.01267350 ; Red_Width=0.071
K16_NIR=0.01542420 ; NIR_Width=0.114
# 1st column: k′ Conversion Factors for 8-Bit products
# Effective Bandwidths [μm] per Spectral Band: as above
K8_Pan10=1.02681367
K8_Pan13=1.02848939
K8_Pan18=1.02794702
K8_Pan24=1.02989685
K8_Pan32=1.02739898
K8_Blue=1.12097834
K8_Green=1.37652632
K8_Red=1.30924587
K8_NIR=0.98368622
# Bands
Spectral_Bands="Pan Blue Green Red NIR"
# loop over all bands
for BAND in ${Spectral_Bands}; do
# which Pan TDI level?
if [[ ${BAND} == Pan ]]
then
eval Pan="Pan${PanTDI}"; echo "Note, processing for ${Pan}"
eval BAND="${Pan}"
fi
# some echo...
echo "Processing the ${BAND} spectral band"
# set band parameters as variables -- using K16 for 16-bit data!
eval K_BAND="K16_${BAND}"
eval BAND_Width="${BAND}_Width"
echo "Band Parameters set to K=${!K_BAND}, Bandwidth=${!BAND_Width}"
if [[ ${BAND} == Pan* ]]
then
BAND="Pan" ; echo "Processing the ${BAND} spectral band"
fi
# set region
g.region rast=${BAND}_DNs #-pg
echo "Region matching the ${BAND} spectral band"
# conversion to Radiance based on (1) -- attention: 32-bit calculations required
r.mapcalc "${BAND}_Radiance = ( double(${!K_BAND}) * ${BAND}_DNs ) / ${!BAND_Width}"
r.info -r "${BAND}_Radiance"
# add info
r.support map=${BAND}_Radiance \
title="" \
units="W / sq.m. / μm / ster" \
description="Top-of-Atmosphere `echo ${BAND}` band spectral Radiance [W/m^2/sr/μm]" \
source1='"Radiometric Use of QuickBird Imagery, Technical Note (2005)," by Keith Krause, Digital Globe'
eval BAND_Esun="${BAND}_Esun"; echo "Using Esun=${!BAND_Esun}"
# calculate ToAR -- ${BAND}_Radiance is already 32-bit -- see above!
r.mapcalc "${BAND}_ToAR = \
( ${PI} * ${BAND}_Radiance * ${ESD}^2 ) / ( ${!BAND_Esun} * cos(${SZA}) )"
r.info -r ${BAND}_ToAR
# add some metadata
r.support map=${BAND}_ToAR \
title="echo ${BAND} band (Top of Atmosphere Reflectance)" \
units="Unitless planetary reflectance" \
description="Top of Atmosphere `echo ${BAND}` band spectral Reflectance (unitless)" \
source1='"Radiometric Use of QuickBird Imagery, Technical Note (2005)," by Keith Krause, Digital Globe' \
source2="Digital Globe" \
history="PI=3.14159265358; K=${!K_BAND}; Bandwidth=${!BAND_Width}; ESD=${ESD}; Esun=${!BAND_Esun}; SZA=${SZA}"
done
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.
Pan Sharpening
In GRASS ver. 6.x, one can use the i.fusion.brovey module.
i.fusion.brovey -q ms1=qbird_green ms2=qbird_nir ms3=qbird_red pan=qbird_pan outputprefix=brov