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 ${\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

Modules overview

Outside of GRASS

• gdalinfo
• gdalwarp

Inside GRASS GIS

• r.in.gdal
• r.mapcalc
• r.colors
• i.landsat.rgb
• i.pansharpen
• i.vi
• i.segment
• ...

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:

Followig needs cleaning

# 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}) )"

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

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