[GRASS-user] Filtering high "outliers" in Landsat reflectance imagery?

NikosAlexandris · February 20, 2013, 11:54am

Advanced users,

may I seek for some recommendation on filtering Landsat reflectance outliers?

I have many Landsat scenes pre-processed (DN to Radiance/Reflectance, Cloud-
masked, Topo-Corrected, band-wise patched in large maps over Greece) and ready
for further explorations.

Applying "color=grey -e OR color=grey1.0 -e" doesn't work well for all larger
maps (after patching -- call it mosaicking if you prefer). That is, some maps
appear too dark (i.e. band 3).

I guess that this may be due to abnormally (?) high reflectance values. I
guess those are artefacts, or not?

The univariate stats of 1+6 bands look fine to me, e.g.:

mean: 298.591
mean of absolute values: 298.591 ### this is Temperature in K

mean: 0.0453416
mean of absolute values: 0.0453416

mean: 0.0490654
mean of absolute values: 0.0490654

mean: 0.0785879
mean of absolute values: 0.0785879

mean: 0.139015
mean of absolute values: 0.139015

mean: 0.121867
mean of absolute values: 0.121867

mean: 0.0845493
mean of absolute values: 0.0845493

Yet, the max Top-of-Canopy Reflectances:

maximum: 326.271 # This is Temperature in K
maximum: 172.05
maximum: 117.96
maximum: 775.934
maximum: 1.66005
maximum: 120.506
maximum: 477.744

Is my understanding correct?
Is there a safe criterion to filter high reflectance values? Could they be
attributed to other sources, e.g. fires?

Can I use some different color rules/scheme which will "ignore" too high
reflectances? Simply "color=grey1.0" or other based on stddev, quantiles?

Thank you in advance for your invaluable time,
Nikos

NikosAlexandris · February 20, 2013, 2:35pm

On Wednesday 20 of February 2013 13:54:11 Nikos Alexandris wrote:

Advanced users,

may I seek for some recommendation on filtering Landsat reflectance
outliers?

(Thanks Yann)

Reflectance is a unit-less ratio, ranging from 0 to 1.0. So I can flatten
everything >1.0 to 1.0.

Best, Nikos

neteler · February 23, 2013, 10:35am

On Wed, Feb 20, 2013 at 3:35 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

On Wednesday 20 of February 2013 13:54:11 Nikos Alexandris wrote:

Advanced users,

may I seek for some recommendation on filtering Landsat reflectance
outliers?

(Thanks Yann)

Reflectance is a unit-less ratio, ranging from 0 to 1.0. So I can flatten
everything >1.0 to 1.0.

I would run r.neighbors with a 3x3 or 5x5 moving window and averaging
or likewise. The question is always: are they outliers or data?

Markus

NikosAlexandris · February 23, 2013, 12:17pm

Nikos A wrote:

>> Advanced users, may I seek for some recommendation on filtering Landsat
>> reflectance outliers?

Nikos A (replying to self):

> (Thanks Yann)
> Reflectance is a unit-less ratio, ranging from 0 to 1.0. So I can flatten
> everything >1.0 to 1.0.

Markus N:

I would run r.neighbors with a 3x3 or 5x5 moving window and averaging
or likewise. The question is always: are they outliers or data?

Wild guessing: too high to be data!?

Nikos

NikosAlexandris · July 29, 2013, 10:55pm

Nikos A wrote:

> >> Advanced users, may I seek for some recommendation on filtering Landsat
> >> reflectance outliers?

Nikos A (replying to self):

> > (Thanks Yann)
> > Reflectance is a unit-less ratio, ranging from 0 to 1.0. So I can
> > flatten everything >1.0 to 1.0.

Markus N:

> I would run r.neighbors with a 3x3 or 5x5 moving window and averaging
> or likewise. The question is always: are they outliers or data?

Wild guessing: too high to be data!?
Nikos

Coming back to this! I think that Markus' suggestion was/is "correct".

However, at the time of having to decide what to do, I simply flattened them.
They were only a few though, not a real issue.

Anyhow, just for completeness an "official" definition of reflectance, here
is one copy-pasted from
<http://landsat.gsfc.nasa.gov/references/glossary.html>:

--%<---
reflectance - A measure of the ratio of outgoing to incoming radiation
calculated by converting a radiometrically calibrated image to an innate
characteristic of the target being observed. Calibrated at-satellite spectral
radiance is converted to unitless reflectance by separating out the
atmospheric component of the reflective band radiance and assuming that the
target is a Lambertian reflector, re-radiating incident solar radiation
equally in all directions. In general, reflectance is a function of incident
angle of the energy, viewing angle of the sensor, spectral wavelength, and
bandwidth, and the nature of the object. Also see planetary albedo,
bidirectional reflectance and atmospheric correction. (Source: Dr. John
Barker)
--->%--

Nikos

Eric_Goddard · August 1, 2013, 3:13pm

Hi, I need to do the same thing to some WV-2 imagery. I'm using the command

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)",
output=wv2_out, input_rast=wv2_in) from python (input and output names
changed for simplicity) but instead of replacing the values greater
than 1 with 1, it replaces them with NaNs. Is there something wrong
with my mapcalc expession?

On Mon, Jul 29, 2013 at 5:55 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

Nikos A wrote:

> >> Advanced users, may I seek for some recommendation on filtering Landsat
> >> reflectance outliers?

Nikos A (replying to self):

> > (Thanks Yann)
> > Reflectance is a unit-less ratio, ranging from 0 to 1.0. So I can
> > flatten everything >1.0 to 1.0.

Markus N:

> I would run r.neighbors with a 3x3 or 5x5 moving window and averaging
> or likewise. The question is always: are they outliers or data?

Wild guessing: too high to be data!?
Nikos

Coming back to this! I think that Markus' suggestion was/is "correct".

However, at the time of having to decide what to do, I simply flattened them.
They were only a few though, not a real issue.

Anyhow, just for completeness an "official" definition of reflectance, here
is one copy-pasted from
<http://landsat.gsfc.nasa.gov/references/glossary.html>:

--%<---
reflectance - A measure of the ratio of outgoing to incoming radiation
calculated by converting a radiometrically calibrated image to an innate
characteristic of the target being observed. Calibrated at-satellite spectral
radiance is converted to unitless reflectance by separating out the
atmospheric component of the reflective band radiance and assuming that the
target is a Lambertian reflector, re-radiating incident solar radiation
equally in all directions. In general, reflectance is a function of incident
angle of the energy, viewing angle of the sensor, spectral wavelength, and
bandwidth, and the nature of the object. Also see planetary albedo,
bidirectional reflectance and atmospheric correction. (Source: Dr. John
Barker)
--->%--

Nikos
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

NikosAlexandris · August 1, 2013, 8:07pm

Eric Goddard wrote:

Hi, I need to do the same thing to some WV-2 imagery. I'm using the command

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)",
output=wv2_out, input_rast=wv2_in)

The expression seems fine -- in general. I also see in my *bash* history
entries like

--%<---
r.mapcalc "TopoCorr.B.Trimmed.DOS1.ToCR.1 = if( TopoCorr.B.Trimmed.DOS1.ToCR.1

1.0, 1.0, TopoCorr.B.Trimmed.DOS1.ToCR.1 )"

or

for MAP in `g.mlist rast pat=TopoCorr.B.Trimmed.DOS1.ToCR.?`; do r.mapcalc
"${MAP} = if( ${MAP} > 1.0, 1.0, ${MAP} )"; done
--->%--

Back to yours:

# in python, you instruct
grass.mapcalc(

# and put code inside
"$output = if($input_rast>1.0, 1.0, $input_rast

# then close it
)"

Maybe it should be (not tested :-p):

grass.mapcalc( "$output = if( $input_rast > 1.0, 1.0, $input_rast)" )

?

Nikos

[rest deleted]

NikosAlexandris · August 1, 2013, 8:19pm

On Thursday 01 of August 2013 10:13:56 Eric Goddard wrote:

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)",
output=wv2_out, input_rast=wv2_in)

Hmm.. sorry, the expression is correct. Seems odd. Is the $input_rast really
"there"?

Nikos

Eric_Goddard · August 1, 2013, 8:48pm

Yes, its there. All the values that are less than 1 are correct from
the $input_rast. Strangely enough, I just tested it on the same image
on my home machine and it worked correctly. Both machines running
grass7.

On Thu, Aug 1, 2013 at 3:19 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

On Thursday 01 of August 2013 10:13:56 Eric Goddard wrote:

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)",
output=wv2_out, input_rast=wv2_in)

Hmm.. sorry, the expression is correct. Seems odd. Is the $input_rast really
"there"?

Nikos

Glynn_Clements1 · August 2, 2013, 2:21am

Eric Goddard wrote:

Hi, I need to do the same thing to some WV-2 imagery. I'm using the command

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)", output=wv2_out, input_rast=wv2_in)

from python (input and output names
changed for simplicity) but instead of replacing the values greater
than 1 with 1, it replaces them with NaNs. Is there something wrong
with my mapcalc expession?

No.

Is there anything unusual about the actual map names? Names containing
any of the characters

^ # ( ) + - % > < ! & | ? : ; ~

need to quoted.

--
Glynn Clements <glynn@gclements.plus.com>

Eric_Goddard · August 2, 2013, 1:38pm

On Thu, Aug 1, 2013 at 9:21 PM, Glynn Clements <glynn@gclements.plus.com> wrote:

Eric Goddard wrote:

Hi, I need to do the same thing to some WV-2 imagery. I'm using the command

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)", output=wv2_out, input_rast=wv2_in)

from python (input and output names
changed for simplicity) but instead of replacing the values greater
than 1 with 1, it replaces them with NaNs. Is there something wrong
with my mapcalc expession?

No.

Is there anything unusual about the actual map names? Names containing
any of the characters

^ # ( ) + - % > < ! & | ? : ; ~

need to quoted.

No special characters. An example image name is
x2NOV04172233_M3DS_R6C2_052823926030_01_P001_TOAR.1. The mapcalc
evaluates successfully (and correctly!) on my machine at home. My work
machine is using the following build:

version=7.0.svn
date=2013
revision=56943
build_date=2013-08-01

./configure --prefix=/usr/local/grass70 --with-postgres=yes
--with-postgres-includes=/usr/include/postgresql --enable-largefile
--with-sqlite --with-freetype=yes
--with-freetype-includes=/usr/include/freetype2
--with-proj-share=/usr/share/proj --with-python --with-cxx
--with-wxwidgets --with-tcltk-includes=/usr/include/tcl8.5
--with-odbc=yes --enable-64bit --with-geos=yes --with-openmp=yes
--with-liblas=yes --with-cairo=yes --with-opencl=yes --with-pthread

I'll have to bost my home machine's version later today. I compiled
the latest svn yesterday to see if that would resolve the issue, but I
still get null values. Not sure which version I was using before. Any
recommendations on how I may be able to track the error down? I'm
using Ubuntu 13.04 64-bit. I wonder if my issue could be linked to
http://osgeo-org.1560.x6.nabble.com/seemingly-random-errors-in-r-mapcalc-td5070327.html
? Unfortunately my mapcalc operation doesn't actually raise an error
so I don't think Soren's recommendation would help in my case.

Eric_Goddard · August 2, 2013, 7:50pm

OK, I discovered what was different between my work and home machines.
The mapcalc if statement doesn't evaluate properly if using
r.external.out. If r.external isn't used, it evaluates correctly. Bug?

On Fri, Aug 2, 2013 at 8:38 AM, Eric Goddard <egoddard1010@gmail.com> wrote:

On Thu, Aug 1, 2013 at 9:21 PM, Glynn Clements <glynn@gclements.plus.com> wrote:

Eric Goddard wrote:

Hi, I need to do the same thing to some WV-2 imagery. I'm using the command

grass.mapcalc("$output = if($input_rast>1.0, 1.0, $input_rast)", output=wv2_out, input_rast=wv2_in)

from python (input and output names
changed for simplicity) but instead of replacing the values greater
than 1 with 1, it replaces them with NaNs. Is there something wrong
with my mapcalc expession?

No.

Is there anything unusual about the actual map names? Names containing
any of the characters

^ # ( ) + - % > < ! & | ? : ; ~

need to quoted.

No special characters. An example image name is
x2NOV04172233_M3DS_R6C2_052823926030_01_P001_TOAR.1. The mapcalc
evaluates successfully (and correctly!) on my machine at home. My work
machine is using the following build:

version=7.0.svn
date=2013
revision=56943
build_date=2013-08-01

./configure --prefix=/usr/local/grass70 --with-postgres=yes
--with-postgres-includes=/usr/include/postgresql --enable-largefile
--with-sqlite --with-freetype=yes
--with-freetype-includes=/usr/include/freetype2
--with-proj-share=/usr/share/proj --with-python --with-cxx
--with-wxwidgets --with-tcltk-includes=/usr/include/tcl8.5
--with-odbc=yes --enable-64bit --with-geos=yes --with-openmp=yes
--with-liblas=yes --with-cairo=yes --with-opencl=yes --with-pthread

I'll have to bost my home machine's version later today. I compiled
the latest svn yesterday to see if that would resolve the issue, but I
still get null values. Not sure which version I was using before. Any
recommendations on how I may be able to track the error down? I'm
using Ubuntu 13.04 64-bit. I wonder if my issue could be linked to
http://osgeo-org.1560.x6.nabble.com/seemingly-random-errors-in-r-mapcalc-td5070327.html
? Unfortunately my mapcalc operation doesn't actually raise an error
so I don't think Soren's recommendation would help in my case.

Glynn_Clements1 · August 3, 2013, 3:22pm

Eric Goddard wrote:

OK, I discovered what was different between my work and home machines.
The mapcalc if statement doesn't evaluate properly if using
r.external.out. If r.external isn't used, it evaluates correctly. Bug?

What parameters are you using for r.external.out?

--
Glynn Clements <glynn@gclements.plus.com>

Eric_Goddard · August 4, 2013, 12:59am

I tried with tif and img files, same results. I didn’t specify any creation options.

r.external.out directory=/data/GIS/satellite_imagery/shelby_wv2/toar extension=.TIF format=GTiff

r.external.out directory=/data/GIS/satellite_imagery/shelby_wv2/toar extension=.IMG format=HFA

On Aug 3, 2013 10:22 AM, “Glynn Clements” <glynn@gclements.plus.com> wrote:

Eric Goddard wrote:

OK, I discovered what was different between my work and home machines.
The mapcalc if statement doesn’t evaluate properly if using
r.external.out. If r.external isn’t used, it evaluates correctly. Bug?

What parameters are you using for r.external.out?

–
Glynn Clements <glynn@gclements.plus.com>

neteler · August 4, 2013, 8:20am

On Sun, Aug 4, 2013 at 2:59 AM, Eric Goddard <egoddard1010@gmail.com> wrote:

I tried with tif and img files, same results. I didn't specify any creation
options.

r.external.out directory=/data/GIS/satellite_imagery/shelby_wv2/toar
extension=.TIF format=GTiff

r.external.out directory=/data/GIS/satellite_imagery/shelby_wv2/toar
extension=.IMG format=HFA

I am not sure that the dot should be given for the "extension" parameter.

Markus

Glynn_Clements1 · August 5, 2013, 1:45am

Eric Goddard wrote:

OK, I discovered what was different between my work and home machines.
The mapcalc if statement doesn't evaluate properly if using
r.external.out. If r.external isn't used, it evaluates correctly. Bug?

Yes. A bug was introduced in r38127 which (by sheer coincidence) only
affects floating-point GDAL output (and primarily affects
double-precision GDAL output). For FCELL, zeros will be mapped to
null; for DCELL, any value where the bottom 32 bits are zero (which
includes all small integers) will be mapped to null.

This should be fixed by r57414.

--
Glynn Clements <glynn@gclements.plus.com>

Eric_Goddard · August 5, 2013, 1:49am

Thanks, looking forward to trying it out in the new version.

On Sun, Aug 4, 2013 at 8:45 PM, Glynn Clements <glynn@gclements.plus.com> wrote:

Eric Goddard wrote:

OK, I discovered what was different between my work and home machines.
The mapcalc if statement doesn't evaluate properly if using
r.external.out. If r.external isn't used, it evaluates correctly. Bug?

Yes. A bug was introduced in r38127 which (by sheer coincidence) only
affects floating-point GDAL output (and primarily affects
double-precision GDAL output). For FCELL, zeros will be mapped to
null; for DCELL, any value where the bottom 32 bits are zero (which
includes all small integers) will be mapped to null.

This should be fixed by r57414.

--
Glynn Clements <glynn@gclements.plus.com>