[GRASS-user] Appropriate input/output range for i.atcorr WorldView2?

Can anyone speak to the appropriate range/datatype of WorldView2 rasters when using the i.atcorr module?

I.atcorr requires the input raster be converted to at-sensor-radiance first, which I have done according to DigitalGlobe documention. DG suggests all transformations be done in FLOAT32 datatype, and i.atcorr requires a specified input/output range that is positive.

I convert from the original DN’s (delivered in UINT16) to radiance (FLOAT32), and while my radiance values are positive, the datatype range is not. In another sense, WV2 is encoded with an 11-bit dynamic range (i.e. 2^11-bit range), but I am unsure if any of these values are right for the module… For example, Sentinel 2 L1C are also encoded as UINT16 with 15 significant bits. However, documentation states that the correct i.atcorr Sentinel range is 1-10000; 10000 being a specified QUANTIFICATION_VALUE by Sentinel.

Is there a WorldView 2 equivalent to this quantification value? Is there an alternative range specific to WV2? Might this range be affected by the radiance conversion? Respectively, what would be the appropriate output/rescale range?

In general- Is the i.atcorr module seeking the literal range of the image pixel values, the range of the datatype they’re in, or the range of data the satellite can collect?

···

Paige Byassee
Appalachian State University '19
B.S., Ecology|Certificate in GIS
(704)488-0872
byasseepaige@gmail.com

Hi Paige,

On Tue, Apr 7, 2020 at 1:45 AM Paige Byassee <byasseepaige@gmail.com> wrote:

Can anyone speak to the appropriate range/datatype of WorldView2 rasters when using the i.atcorr module?

I.atcorr requires the input raster be converted to at-sensor-radiance first, which I have done according to DigitalGlobe documention.
DG suggests all transformations be done in FLOAT32 datatype, and i.atcorr requires a specified input/output range that is positive.

Yes, as per manual
https://grass.osgeo.org/grass78/manuals/i.atcorr.html:

Input range
  range=min,max
  Default: 0,255

.. the default range is 8bit. Since you use WorldView2, the range
parameter needs to be modified (see below).

I have a WorldView2 scene available here, so I am looking into the
details with gdalinfo:

Metadata:
  AREA_OR_POINT=Area
  METADATATYPE=DG
  TIFFTAG_COPYRIGHT=(C) COPYRIGHT 2018 DigitalGlobe, Inc., Longmont CO USA 80503
  TIFFTAG_DATETIME=2018:12:20 09:13:25
  TIFFTAG_IMAGEDESCRIPTION={
  bandList =
  [
    6;
...
    9;
  ]
}
  TIFFTAG_MAXSAMPLEVALUE=2047
  TIFFTAG_MINSAMPLEVALUE=0
Image Structure Metadata:
  INTERLEAVE=BAND
RPC Metadata:
...
Band 1 Block=3217x162 Type=UInt16, ColorInterp=Red
  Overviews: 1609x1660, 805x830, 403x415, 202x208
Band 2 Block=3217x162 Type=UInt16, ColorInterp=Green
...

So, UInt16 also here. Now:

I convert from the original DN's (delivered in UINT16) to radiance (FLOAT32),

I am not entirely convinced that you need to convert UINT16 to FLOAT32
like that.
I believe that you need to apply the formula along with the respective
band calibration values ABSCALFACTOR and EFFECTIVEBANDWIDTH (the same
is also stored in the IMD files) from the XML metadata files:

grep 'ABSCALFACTOR\|EFFECTIVEBANDWIDTH'
058891334020_01_P001_MUL/18DEC03184338-M2AS-058891334020_01_P001.XML
            <ABSCALFACTOR>9.295654000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>4.730000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>9.748051000000001e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>5.430000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>7.541495000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>6.300000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>5.101088000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>3.740000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>1.103623000000000e-02</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>5.740000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>4.539619000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>3.930000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>1.224380000000000e-02</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>9.890000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>9.042234000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>9.959999999999999e-02</EFFECTIVEBANDWIDTH>

Citing from "Absolute Radiometric Calibration, Prepared By: Michele A. Kuester"
https://dg-cms-uploads-production.s3.amazonaws.com/uploads/document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf

  "The top-of-atmosphere radiance, L, in units of Wμm-1m-2sr-1, is
then found from the DigitalGlobe
    image product for each band by converting from digital numbers
(DN) using the equation,

    L = Gain * DN * (abscalfactor/effective bandwidth) + Offset

    The TDI specific "abscalfactor" and "effectiveBandwidth" are
delivered with the imagery in the metadata file.
    The digital number, DN, is the pixel value found in the imagery.
The Gain and Offset are the absolute
    radiometric calibration band dependent adjustment factors that are
given in Table 1. Note that these
    are not necessarily stagnant values and they are revisited annually.
"

(the "Table 1" is found in the same PDF file above. You may want to
check if a newer table version exists).

You can apply this formula with r.mapcalc, for each band to obtain TOA
data from DN.

and while my radiance values are positive, the datatype range is not. In another sense, WV2 is encoded with an 11-bit dynamic range (i.e. 2^11-bit range), but I am unsure if any of these values are right for the module... For example, Sentinel 2 L1C are also encoded as UINT16 with 15 significant bits. However, documentation states that the correct i.atcorr Sentinel range is 1-10000; 10000 being a specified QUANTIFICATION_VALUE by Sentinel.

Is there a WorldView 2 equivalent to this quantification value? Is there an alternative range specific to WV2? Might this range be affected by the radiance conversion? Respectively, what would be the appropriate output/rescale range?

I hope the explanations above help.
Please let us know how it goes.

In general- Is the i.atcorr module seeking the literal range of the image pixel values, t
he range of the datatype they're in, or the range of data the satellite can collect?

Well, unfortunately i.atcorr isn't yet clever enough to parse the
values itself from the metadata. So, one needs to find the values and
apply them.

Here is an overview:
https://grasswiki.osgeo.org/wiki/Atmospheric_correction

and this page should IMHO be updated (your input is most welcome here
- just register there to edit or comment here)
https://grasswiki.osgeo.org/wiki/WorldView#Spectral_Radiance

Best,
Markus

--
Markus Neteler, PhD
https://www.mundialis.de - free data with free software
https://grass.osgeo.org
https://courses.neteler.org/blog

On 2020-04-08 15:17, Markus Neteler wrote:
..

I believe that you need to apply the formula along with the respective
band calibration values ABSCALFACTOR and EFFECTIVEBANDWIDTH (the same
is also stored in the IMD files) from the XML metadata files:

grep 'ABSCALFACTOR\|EFFECTIVEBANDWIDTH'
058891334020_01_P001_MUL/18DEC03184338-M2AS-058891334020_01_P001.XML
            <ABSCALFACTOR>9.295654000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>4.730000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>9.748051000000001e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>5.430000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>7.541495000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>6.300000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>5.101088000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>3.740000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>1.103623000000000e-02</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>5.740000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>4.539619000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>3.930000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>1.224380000000000e-02</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>9.890000000000000e-02</EFFECTIVEBANDWIDTH>
            <ABSCALFACTOR>9.042234000000000e-03</ABSCALFACTOR>
            <EFFECTIVEBANDWIDTH>9.959999999999999e-02</EFFECTIVEBANDWIDTH>

Citing from "Absolute Radiometric Calibration, Prepared By: Michele A. Kuester"
https://dg-cms-uploads-production.s3.amazonaws.com/uploads/document/file/209/ABSRADCAL_FLEET_2016v0_Rel20170606.pdf

  "The top-of-atmosphere radiance, L, in units of Wμm-1m-2sr-1, is
then found from the DigitalGlobe
    image product for each band by converting from digital numbers
(DN) using the equation,

    L = Gain * DN * (abscalfactor/effective bandwidth) + Offset

    The TDI specific "abscalfactor" and "effectiveBandwidth" are
delivered with the imagery in the metadata file.
    The digital number, DN, is the pixel value found in the imagery.
The Gain and Offset are the absolute
    radiometric calibration band dependent adjustment factors that are
given in Table 1. Note that these
    are not necessarily stagnant values and they are revisited annually.
"

(the "Table 1" is found in the same PDF file above. You may want to
check if a newer table version exists).

You can apply this formula with r.mapcalc, for each band to obtain TOA
data from DN.

See also https://gitlab.com/NikosAlexandris/i.worldview.toar or https://github.com/NikosAlexandris/i.worldview.toar.

Maybe useful (though ratehr badly programmed at the time).

Nikos

Thank you Nikos & Marcus,

Converting to TOA radiance (or reflectance) is not really the issue… I have already been following “Absolute Radiometric Calibration, Prepared By: Michele A. Kuester” using the formula for L (done in R, screenshot attached). I do think it is necessary that this be done in FLOAT32 to retain enough detail and per DigitalGlobe recommendations.

So, if I covert to TOA radiance as you’ve described, I am left with 32-bit pixels. I then import the images using r.in.gdal and launch i.atcorr. Would the input range then be 0, 4294967295 (i.e. (2^32) - 1)?

If so, I wonder what a good output range might be, as I presume using r.out.gdal on a 32-bit image could take a long time.

Thank you for your help,
Paige

(attachments)

RadCode.PNG

···

Paige Byassee
Appalachian State University '19
B.S., Ecology|Certificate in GIS
(704)488-0872
byasseepaige@gmail.com

Hello,

I am trying to project data from ll wgas84 location to custom projection. I end up with: ERROR: Unable to open element file <> for <WIND@Morava>. Data were imported to wgs 84 location withou any problems (r.in.gdal) and than I made smaller area using r.resample.

I am using Grass 7.8.2 on Ubuntu 18.04.

Cheers,

Pavel