[GRASS-dev] On (Landsat) imagery naming patterns

Nikos_Alexandris · February 20, 2013, 11:31am

Greetings to the dev-list.

I successfully tested <i.histo.match>, by feeding it with 12 Landsat scenes
(working on identical bands, different Path/Row tiles) in one go. It works and
I really enjoy it.

However, there is some "problem":

<i.histo.match> wont accept the "." character as part of the image/raster map
names. This is an SQL restriction (thanks Luca D for the quick help) as the
module utilises sqlite internally for the calculations.

Although I am familiar with this restriction, due to pressing deadlines, I
overlooked that "point" while designing my workflow and decided to go for a
specific naming pattern for all Landsat scenes at each (pre-)processing step
(working currently on 25+ tiles, for summer and winter that cover the entire
Greek territory).

In addition, the respective manuals for modules like <i.landsat.toar>,
<i.landsat.acca>, <i.topo.corr> hold in their examples/descriptions (thus,
someone can easily interpret that these modules "encourage" the use of) names
like "B.", "B.toar", "toar.5", et.c.

I certainly like to keep, while on intermediate processing steps, long names
whose parts are separated by either the underscore character ( _ ) or the dot
character ( . ).

This is, from my point of view, an inconsistency which adds-up some confusion
-- especially to new users who will try/need to use <i.histo.match> at some
point.

As changing the module or adding internally a solution will add complexity, I
guess it is better to adhere to some kind of "strict" naming rules/patterns in
the manuals' examples, in the wiki and elsewhere?

In time I will try to propose new wording and examples for all Landsat related
modules/manuals I have used/read.

GRASS has (almost) everything required for a "perfect" Landsat processing
work-flow. This small inconsistency breaks the harmony :-).

Best, N

Nikos_Alexandris · February 20, 2013, 2:27pm

On Wednesday 20 of February 2013 13:31:07 Nikos Alexandris wrote:

Greetings to the dev-list.

I successfully tested <i.histo.match>, by feeding it with 12 Landsat scenes
(working on identical bands, different Path/Row tiles) in one go. It works
and I really enjoy it.

Could it be that histo-matching works nicely only on 0-255 ranges but doesn't
play nice in ranges 0-1.0? Integers only?

For example, when feeding images ranging from 0-1.0 (all Landsat bands 3) with
the optional parameter "max=1", many of the histo-matched images are "flatten"
to 1.

I will try to "recode" 0-1.0 images to 0-255 and repeat the histo-matching
process.

Best, Nikos

Markus_Metz · February 20, 2013, 2:54pm

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

On Wednesday 20 of February 2013 13:31:07 Nikos Alexandris wrote:

Greetings to the dev-list.

I successfully tested <i.histo.match>, by feeding it with 12 Landsat scenes
(working on identical bands, different Path/Row tiles) in one go. It works
and I really enjoy it.

Could it be that histo-matching works nicely only on 0-255 ranges but doesn't
play nice in ranges 0-1.0? Integers only?

Yes, integers only. The i.histo.match script forces input raster
values to integers. But it might not be restricted to the 0-255 range.

For example, when feeding images ranging from 0-1.0 (all Landsat bands 3) with
the optional parameter "max=1", many of the histo-matched images are "flatten"
to 1.

I will try to "recode" 0-1.0 images to 0-255 and repeat the histo-matching
process.

For rescaling to 0-1000, you could try something like

r.mapcalc "bandX_int = round(1000 * bandX)"

HTH,

Markus M

Best, Nikos

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

Nikos_Alexandris · February 20, 2013, 3:19pm

Nikos Alexandris:

>> Greetings to the dev-list.
>> I successfully tested <i.histo.match>, by feeding it with 12 Landsat
>> scenes (working on identical bands, different Path/Row tiles) in one go.
>> It works and I really enjoy it.

Nikos A:

> Could it be that histo-matching works nicely only on 0-255 ranges but
> doesn't play nice in ranges 0-1.0? Integers only?

Markus Metz wrote:

Yes, integers only. The i.histo.match script forces input raster
values to integers. But it might not be restricted to the 0-255 range.

(
Oh how nice it would be to do the tranformation internally
)

> For example, when feeding images ranging from 0-1.0 (all Landsat bands 3)
> with the optional parameter "max=1", many of the histo-matched images are
> "flatten" to 1. I will try to "recode" 0-1.0 images to 0-255 and repeat
> the histo-matching process.

For rescaling to 0-1000, you could try something like
r.mapcalc "bandX_int = round(1000 * bandX)"

Nice! But I wonder why round-ing is required? I want to jump back afterwards
(i.e. divide by 1000 the histo-matched images). Rounding does not hurt much,
does it?

Thank you Markus, Nikos

Markus_Metz · February 20, 2013, 3:34pm

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

Nikos Alexandris:

Markus Metz wrote:

For rescaling to 0-1000, you could try something like
r.mapcalc "bandX_int = round(1000 * bandX)"

Nice! But I wonder why round-ing is required?

rounding is not required, just a bit more accurate than truncating
with int(1000 * bandX).

Markus M

Nikos_Alexandris · February 20, 2013, 3:50pm

MM:

>> For rescaling to 0-1000, you could try something like
>> r.mapcalc "bandX_int = round(1000 * bandX)"

NA:

> Nice! But I wonder why round-ing is required?

MM:

rounding is not required, just a bit more accurate than truncating
with int(1000 * bandX).

Great Wizard, thnx

N

ps- note-to-self: add this in Wiki when writing about histogram matching of
Landsat reflectances.

Yann3 · February 20, 2013, 5:53pm

Hmmm g.rename originalnname with newname (sed substitute “.” with “_” in originalname)

I agree this additional step is less elegant in the work flow
Y

On Feb 20, 2013 12:31 PM, “Nikos Alexandris” <nik@nikosalexandris.net> wrote:

Greetings to the dev-list.

I successfully tested <i.histo.match>, by feeding it with 12 Landsat scenes
(working on identical bands, different Path/Row tiles) in one go. It works and
I really enjoy it.

However, there is some “problem”:

<i.histo.match> wont accept the “.” character as part of the image/raster map
names. This is an SQL restriction (thanks Luca D for the quick help) as the
module utilises sqlite internally for the calculations.

Although I am familiar with this restriction, due to pressing deadlines, I
overlooked that “point” while designing my workflow and decided to go for a
specific naming pattern for all Landsat scenes at each (pre-)processing step
(working currently on 25+ tiles, for summer and winter that cover the entire
Greek territory).

In addition, the respective manuals for modules like <i.landsat.toar>,
<i.landsat.acca>, <i.topo.corr> hold in their examples/descriptions (thus,
someone can easily interpret that these modules “encourage” the use of) names
like “B.”, “B.toar”, “toar.5”, et.c.

I certainly like to keep, while on intermediate processing steps, long names
whose parts are separated by either the underscore character ( _ ) or the dot
character ( . ).

This is, from my point of view, an inconsistency which adds-up some confusion
– especially to new users who will try/need to use <i.histo.match> at some
point.

As changing the module or adding internally a solution will add complexity, I
guess it is better to adhere to some kind of “strict” naming rules/patterns in
the manuals’ examples, in the wiki and elsewhere?

In time I will try to propose new wording and examples for all Landsat related
modules/manuals I have used/read.

GRASS has (almost) everything required for a “perfect” Landsat processing
work-flow. This small inconsistency breaks the harmony :-).

Best, N

Nikos_Alexandris · February 20, 2013, 7:39pm

On Wednesday 20 of February 2013 23:23:31 Yann Chemin wrote:

Hmmm g.rename originalnname with newname (sed substitute "." with "_" in
originalname)

I agree this additional step is less elegant in the work flow
Y

Hi Yann.

This is, of course, what I did/do for now. I can't change my naming pattern
now anyway. But I think as a long-term goal, the modules should follow
homogeneous naming conventions, as much as possible.

For the records, I use(d) in a script

tr "." "_"

as I have placed the "." as a separator everywhere. I hope this is not
dangerous.

Best regards, Nikos

Hamish · February 20, 2013, 8:36pm

Nikos wrote:

This is, of course, what I did/do for now. I can't
change my naming pattern
now anyway. But I think as a long-term
goal, the modules should follow
homogeneous naming conventions, as much as possible.

It should be easy to make i.histo.match.py do name cleansing and
then back to the orginal, internally. It could make a table/
dictionary of map names vs. img_1, img_2, ... safe column names
then work with the safe names. Make the program[mer] do the hard
work, not the end user..

best,
Hamish

Nikos_Alexandris · February 20, 2013, 11:02pm

Nikos wrote:

> This is, of course, what I did/do for now. I can't
> change my naming pattern
> now anyway. But I think as a long-term
> goal, the modules should follow
> homogeneous naming conventions, as much as possible.

Hamish wrote:

It should be easy to make i.histo.match.py do name cleansing and
then back to the orginal, internally. It could make a table/
dictionary of map names vs. img_1, img_2, ... safe column names
then work with the safe names. Make the program[mer] do the hard
work, not the end user..

(
note-to-self: maybe that would be (another) good exercise for me (later, on
March)?
)

Note, <i.histo.match> itself adds a user-defined suffix which _always_ begins
with a "."

thnx, Nikos

Nikos_Alexandris · February 21, 2013, 9:18am

MM:

> >> For rescaling to 0-1000, you could try something like
> >> r.mapcalc "bandX_int = round(1000 * bandX)"

NA:

> > Nice! But I wonder why round-ing is required?

MM:

> rounding is not required, just a bit more accurate than truncating
> with int(1000 * bandX).

Just FYI, results look nice! I even convert back to 0-1.0 via

r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"

Thnx, Nikos

ps- I wonder if *10000 is *better* for higher precision?

Hamish · February 22, 2013, 3:17am

Nikos wrote:

Just FYI, results look nice! I even convert back to
0-1.0 via

r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"

...

ps- I wonder if *10000 is *better* for higher precision?

depending on the sensor's 8-bitnesss or not, you can probably calculate by hand now many significant digits are useful. A little bit extra probably doesn't hurt.

1 / 2^8 = 0.00390625

Hamish

Nikos_Alexandris · February 22, 2013, 9:49am

Jorge:

Hi, I quick review the code of i.landsat.* and the programms never use
the dot as a part of the output filenames.
The output names are construct as prefix+suffix not as
prefix+'.'+suffix. Then, the question is more on 'style' than on
programmation. The only change needed are in the HTML when a decission
is taken.

Right!

( I also wrote in the beginning:

Nikos A:
..

In addition, the respective manuals for modules like <i.landsat.toar>,
<i.landsat.acca>, <i.topo.corr> hold in their examples/descriptions (thus,
someone can easily interpret that these modules "encourage" the use of)
names like "B.", "B.toar", "toar.5", et.c.

..
)

It's just that people will likely use what the manuals suggest!

Another quick question: What is the better name to the Landsat 8? tm8 or
ot8 (from OLI/TIRS 8).

(
If my view counts a small little bit: "tm8" rings a bell like "thematic
mapper" which is maybe a bit confusing.
)

I used ot8 in the previous version send to Nikos
but before send to Markus it could be changed.

Sorry Jorge, don't have the time to test it now.

Thank you once again for your efforts Jorge
Best, Nikos

Nikos_Alexandris · February 22, 2013, 12:24pm

Nikos wrote:

> Just FYI, results look nice! I even convert back to 0-1.0 via
> r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"

..

> ps- I wonder if *10000 is *better* for higher precision?

Hamish wrote:

depending on the sensor's 8-bitnesss or not, you can probably calculate by
hand now many significant digits are useful. A little bit extra probably
doesn't hurt.
1 / 2^8 = 0.00390625

- Does this add up (practically a lot?) to the computational cost for "mortal"
machines?

And two more questions on <i.histo.match>:

- is there a reason, besides time, not to make <i.histo.match> understand
floating points?

- the algorithm implemented in <i.histo.match> does take Average pixel
values as reference(s) upon which all computations are based on for the histo-
matched values, right?

--%<--- lines 165 ~ 167 --- --- --- %<---
# calculate new value of pixel_frequency
average = (tot/n_images)
cHist = cHist + int(average)
--->%-- --- --- --- --->%-- --- -- --- --->%---

thnx, Nikos

Nikos_Alexandris · July 31, 2013, 10:38pm

Nikos wrote:

> Just FYI, results look nice! I even convert back to 0-1.0 via

> r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"
> ps- I wonder if *10000 is *better* for higher precision?

By the way, that was Landsat5!

depending on the sensor's 8-bitnesss or not, you can probably calculate by
hand now many significant digits are useful. A little bit extra probably
doesn't hurt.

1 / 2^8 = 0.00390625

I didn't pay too much attention back then as I was under enormous time
pressure. But, I think I didn't miss much of the precision with respect to
the final product's scope. Except if the last 4 (or 3, can't remember) digits
make up a great deal when histo-matching.

Back to present. QuickBird is an 11-bit sensor [1] and data are delivered as
either 8-bit or 16-bit. In this case we have

1 / 2^11 = 0.0004882812

Does that say anything about the significant digits? I want to rescale
QuickBird bands in order to use'em with i.pansharpen. I am thinking to go
from Reflectance (double precision) to integer (as above, only this time I'd
multiply with a number as big as it takes to keep all decimals) and then
rescale to [0, 255].

Makes non-/sense? Any idea how many decimals should be preserved for
analysis? I might be hunting "fine digits" for nothing...

In any case, conversions from DNs to Reflectance should be done in a 32-bit
level [2] (that corresponds to: 1 / 2.328306e-10). The same is done with
Landsat imagery though both L5 and L7 data are delivered as 8-bit [3].

Nikos

---
[1] also mentioned in the document
<http://www.digitalglobe.com/downloads/QuickBird_technote_raduse_v1.pdf>, page
7

[2] same document, page 8: "...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."

[3] http://landsat.usgs.gov/how_is_radiance_calculated.php