QuickBird

Overview

QuickBird is a commercial earth observation satellite. Details about the sensor are provided at Digital Globe's QuickBird Data Sheet. Wikipedia's article on QuickBird provides an informative overview as well.

Quick Facts

• Panchromatic (black and white) and Multi-Spectral imagery
• Multi-spectral bands are ordered as

band 1 = blue
     2 = green
     3 = red
     4 = Near IR

• 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.

• Available pixel resolutions

Product	Black & White	Multispectral	Color & Pan-sharpened
Basic	61 cm to 72 cm as collected	2.44 m to 2.88 m as collected	not available
Standard	60 cm or 70 cm	2.4 m or 2.8 m	60 cm or 70 cm
Orthorectified	60 cm or 70 cm	2.4 m or 2.8 m	60 cm or 70 cm

Notes

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.
• Some sample QuickBird2 data are available at Digital Globe's collection of Imagery Product Samples
• The Global Land Cover Facility (GLCF) provides some QuickBird acqusitions for education or research purposes only.

Modules overview

Outside of GRASS

• gdalinfo
• gdalwarp

Inside GRASS GIS

• r.in.gdal
• r.mapcalc
• r.colors
• i.colors.enhance
• i.fusion.hpf (src)
• i.pansharpen
• i.quickbird.toar
• i.vi
• i.segment

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 ${\frac {W}{m^{2}*sr*nm}}$, 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:
$L_{\lambda {\text{Pixel, Band}}}={\frac {K_{\text{Band}}*q_{\text{Pixel, Band}}}{\Delta \lambda _{\text{Band}}}}$
where:
• $L_{\lambda {\text{Pixel,Band}}}$: Top-of-Atmosphere Spectral Radiance image pixels [W-m-2-sr-1-μm-1]
• $K_{\text{Band}}$: the absolute radiometric calibration factor [W-m-2-sr-1-count-1] for a given band
• $q_{\text{Pixel,Band}}$: radiometrically corrected image pixels [counts or Digital Numbers]
• $\Delta _{\lambda _{\text{Band}}}$: 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:

$\rho _{p}={\frac {\pi *L\lambda *d^{2}}{ESUN\lambda *cos(\Theta _{S})}}$


where:

• $\rho$ - Unitless Planetary Reflectance
• $\pi$ - mathematical constant (3.14159265358)
• $L\lambda$ spectral Radiance at the sensor's aperture, from equation... ToADD
• $d$ - Earth-Sun distance in astronomical units, interpolated values
• $Esun$ - Mean solar exoatmospheric irradiance(s) (W/m2/μm), interpolated values
• $cos(\theta _{s})$ - Solar_zenith_angle from the image acquisition's metadata

Importing data

File Formats & Metadata

• GeoTIFF, NITF 2.1, NITF 2.0
• Image bits / pixel: 8 or 16 bits

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

Importing

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

It is very likely that the Red and Blue bands are swapped (in older datasets), 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 Radiometric Use of QuickBird Imagery, Technical Note. 2005-11-07, by Keith Krause.

# 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 ${\frac {W}{m^{2}*\mu m}}$. 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}) )"

Automating Conversions

Scripting the conversion process avoids repeating the same steps for each band separately.

in Python

A custom python script, performing the operations of interest, might be like i.quickbird.toar (for GRASS 7.x)

in bash

Am example bash shell script, might look like the one that follows. Note, however, 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

See 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

HPFA based sharpening

The following screenshots exemplify the High Pass Filtering Addition fusion technique applied in a fragment of the QuickBird acquisition "04APR05050541-X2AS_R1C1-000000186011_01_P001-Sri_Lanka-Kokilai_Lagoon" which is publicly available via the GLCF. An implementation of this technique for GRASS-GIS is available as a grass-addon i.fusion.hpf (src).

Fragment from the high resolution (0.7m) panchromatic image

Fragment from the low resolution (2.8m) multi-spectral RGB composite image

Fragment from the HPFA Pan-Sharpened RGB composite image. Sharpening performed with default parameters.

Fragment from another HPFA Pan-Sharpened RGB composite image. Sharpening performed with center=low and modulator=min parameters.