[GRASS-dev] Landsat 5TM pre-processing - histogram matching - problem

Hello,

I’m totally new user of GRASS and of course I have some problems. I want to do the pre-processing of Landsat (5TM) images (for further band composite, classification and NDVI) and I was using the tips from http://grasswiki.osgeo.org/wiki/LANDSAT#Pre-Processing I did the i.landsat.toar I wanted to do the i.atcorr as well but I have no idea how to do this, the manual was not helpful for me, so I gave up. I did not do i.topo.corr because the land I’m working on is flat.

And I really want to do “radiometrically normalise → one approach via i.histo.match (in grass 7), also known as relative radiometric normalisation – one approach is the histogram matching technique of two or more raster maps”. I have fragments of original landsat images imported to GRASS and i.histo.match works with them without any problems. But there is a problem when I want to do this on images which I got after i.landsat.toar. The output images are monochromatic-one colour. I’ve used command i.histo.match input=_toar2@konfa,86_toar2@konfa (the input files are those which I got after i.landsat.toar) and then r.colors map=86_toar2.match@konfa color=grey. I thought that those steps from grasswiki should be done “step by step”, so first i.landsat.toar, and then on the output files i.histo.match.

Does anyone know how to do this in a proper/correct way?

Best wishes

Joan

On 04/05/15 11:28, joanna mardas wrote:

Hello,

I'm totally new user of GRASS and of course I have some problems. I want
to do the pre-processing of Landsat (5TM) images (for further band
composite, classification and NDVI) and I was using the tips from
http://grasswiki.osgeo.org/wiki/LANDSAT#Pre-Processing I did the
*i.landsat.toar* I wanted to do the i.atcorr as well but I have no idea
how to do this, the manual was not helpful for me, so I gave up. I did
not do i.topo.corr because the land I'm working on is flat.

And I really want to do "/radiometrically normalise → one approach
via*i.histo.match* (in grass 7), also known as relative radiometric
normalisation -- one approach is the histogram matching technique of two
or more raster maps/". I have fragments of original landsat images
imported to GRASS and i.histo.match works with them without any
problems. But there is a problem when I want to do this on images which
I got after i.landsat.toar. The output images are monochromatic-one
colour. I've used command i.histo.match
input=_toar2@konfa,86_toar2@konfa (the input files are those which I got
after i.landsat.toar) and then r.colors map=86_toar2.match@konfa
color=grey. I thought that those steps from grasswiki should be done
"step by step", so first i.landsat.toar, and then on the output files
i.histo.match.

Does anyone know how to do this in a proper/correct way?

ISTR that i.histo.match works with integer imagerie, i.e. recoded to 0-255. Recode your toar images to that range with r.recode and try to run the result through i.histo.match to see if that helps.

Moritz

joanna mardas wrote:

> Hello,

Hi Joanna,

> I'm totally new user of GRASS and of course I have some problems. I want
> to do the pre-processing of Landsat (5TM) images (for further band
> composite, classification and NDVI) and I was using the tips from
> http://grasswiki.osgeo.org/wiki/LANDSAT#Pre-Processing I did the
> *i.landsat.toar* I wanted to do the i.atcorr as well but I have no idea
> how to do this, the manual was not helpful for me, so I gave up. I did
> not do i.topo.corr because the land I'm working on is flat.

I need some time to fix a python script which does that automatically
(see: <https://github.com/NikosAlexandris/i.landsat.atcorr&gt;, currently
broken :-().

> And I really want to do "/radiometrically normalise → one approach
> via*i.histo.match* (in grass 7), also known as relative radiometric
> normalisation -- one approach is the histogram matching technique of two
> or more raster maps/".

Another, simple approach, is described in:

"Radiometric Use of QuickBird Imagery, Technical Note." 2005-11-07, by
Keith Krause.

It is in section 6 in the above mentioned document. I have a
half-written script for this... :frowning:

> I have fragments of original landsat images
> imported to GRASS and i.histo.match works with them without any
> problems.

That is so because i.histo.match, as Moritz notes below, support for
8-bit raster data.

> But there is a problem when I want to do this on images which
> I got after i.landsat.toar. The output images are monochromatic-one
> colour.

Each band is one single image. No matter its bitness, it is "normal" to
appear "mono-chromatic" -- its "normal" to apply a grey scale.

> I've used command i.histo.match
> input=_toar2@konfa,86_toar2@konfa (the input files are those which I got
> after i.landsat.toar) and then r.colors map=86_toar2.match@konfa
> color=grey. I thought that those steps from grasswiki should be done
> "step by step", so first i.landsat.toar, and then on the output files
> i.histo.match.
>
> Does anyone know how to do this in a proper/correct way?

I wrote this steps some time ago. And I still insist that the correct
way is to not use the Digital Numbers directly. Rather, convert to
reflectance, then do whatever you have to.

Moritz Lennert:

ISTR that i.histo.match works with integer imagerie, i.e. recoded to
0-255. Recode your toar images to that range with r.recode and try to
run the result through i.histo.match to see if that helps.

What I have done in the past,

# if Reflectance > 1, set to 1000, else round and multiply by 1000
r.mapcalc --o "TopoCorr_B_Trimmed_DOS1_ToCR_${Band_No}_${SCENE} = if( TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} > 1, 1000, round( 1000 * TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} ) )"

# histogram matching
...

# convert histo-matched images back to floating point: rescale to 0-1.0
r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"

This came out after Markus Metz' advice, if I recall correctly. I would
add another advice (this one was from Yann Chemin): filter out water. I
think too, now, that i.histo.match will perform better.

Nikos

joanna mardas wrote:

> > I have fragments of original landsat images
> > imported to GRASS and i.histo.match works with them without any
> > problems.

Nikos Alexandris :

That is so because i.histo.match, as Moritz notes below, support for
8-bit raster data.

Correction: support is there for integers, I think. So, not limited to
8-bit only.

Nikos

Hello!

You both are my heros! Thank you Moritz and Nikos!

r.recode worked :smiley: I did i.histo.match on two test images and it looks fine. So now the rest of bands, and then band composite and supervised classification (I’ve found nice tutorial on youtube, so I’m gonna follow it, or I will try semi-automatic classification plugin for qgis).

Nikos, you’ve wrote “Each band is one single image. No matter its bitness, it is “normal” to appear “mono-chromatic” – its “normal” to apply a grey scale.” I know that grey scale is normal, but I didn’t have, let’s say “50 shades of grey” :wink: but only one shade, one color-no scale of colors. Forunately, Moritz was right, thanks a lot! Now it looks normal.

And you wrote also that “the correct way is to not use the Digital Numbers directly. Rather, convert to reflectance, then do whatever you have to.” So application of i.landsat.toar is correct, right?

Sorry for all stupid/basic questions but this is my first serious time with grass and landsat image processing.

Greetings

Joanna

Dnia 4-05-2015 o godz. 21:14 Nikos Alexandris napisa³(a):

joanna mardas wrote:

Hello,

Hi Joanna,

I’m totally new user of GRASS and of course I have some problems. I want

to do the pre-processing of Landsat (5TM) images (for further band

composite, classification and NDVI) and I was using the tips from

http://grasswiki.osgeo.org/wiki/LANDSAT#Pre-Processing I did the

i.landsat.toar I wanted to do the i.atcorr as well but I have no idea

how to do this, the manual was not helpful for me, so I gave up. I did

not do i.topo.corr because the land I’m working on is flat.

I need some time to fix a python script which does that automatically

(see: https://github.com/NikosAlexandris/i.landsat.atcorr, currently

broken :-().

And I really want to do "/radiometrically normalise → one approach

viai.histo.match (in grass 7), also known as relative radiometric

normalisation – one approach is the histogram matching technique of two

or more raster maps/".

Another, simple approach, is described in:

“Radiometric Use of QuickBird Imagery, Technical Note.” 2005-11-07, by

Keith Krause.

It is in section 6 in the above mentioned document. I have a

half-written script for this… :frowning:

I have fragments of original landsat images

imported to GRASS and i.histo.match works with them without any

problems.

That is so because i.histo.match, as Moritz notes below, support for

8-bit raster data.

But there is a problem when I want to do this on images which

I got after i.landsat.toar. The output images are monochromatic-one

colour.

Each band is one single image. No matter its bitness, it is “normal” to

appear “mono-chromatic” – its “normal” to apply a grey scale.

I’ve used command i.histo.match

input=_toar2@konfa,86_toar2@konfa (the input files are those which I got

after i.landsat.toar) and then r.colors map=86_toar2.match@konfa

color=grey. I thought that those steps from grasswiki should be done

“step by step”, so first i.landsat.toar, and then on the output files

i.histo.match.

Does anyone know how to do this in a proper/correct way?

I wrote this steps some time ago. And I still insist that the correct

way is to not use the Digital Numbers directly. Rather, convert to

reflectance, then do whatever you have to.

Moritz Lennert:

ISTR that i.histo.match works with integer imagerie, i.e. recoded to

0-255. Recode your toar images to that range with r.recode and try to

run the result through i.histo.match to see if that helps.

What I have done in the past,

if Reflectance > 1, set to 1000, else round and multiply by 1000

r.mapcalc --o "TopoCorr_B_Trimmed_DOS1_ToCR_${Band_No}_${SCENE} = if(

TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} > 1, 1000, round( 1000

  • TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} ) )"

histogram matching

convert histo-matched images back to floating point: rescale to 0-1.0

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

This came out after Markus Metz’ advice, if I recall correctly. I would

add another advice (this one was from Yann Chemin): filter out water. I

think too, now, that i.histo.match will perform better.

Nikos

* joanna mardas <joanna.mardas@wp.pl> [2015-05-05 00:19:02 +0200]:

   Hello!

   You both are my heros! Thank you Moritz and Nikos!

   r.recode worked :smiley: I did i.histo.match on two test images and it looks
   fine.

If you go all the way from DN to Top-of-Atmosphere reflectances (via
i.landsat.toar), then to Top-of-Canopy (via i.atcorr), you'll have
floating point values ranging in [0, 1.0]. If you recode this back to 8-bit,
you should consider whether there is an important loss of the radiometric resolution.

So, recoding to another range is different than converting to integer
numbers and trying to preserve the range.

   So now the rest of bands, and then band composite and supervised
   classification (I've found nice tutorial on youtube, so I'm gonna follow
   it, or I will try semi-automatic classification plugin for qgis).

   Nikos, you've wrote "Each band is one single image. No matter its bitness,
   it is "normal" to appear "mono-chromatic" -- its "normal" to apply a grey
   scale." I know that grey scale is normal, but I didn't have, let's say "50
   shades of grey" :wink: but only one shade, one color-no scale of colors.

Thanks for the details Joanna. I now understand better. After each
processing you can check the range of your image(s) via `r.info -r` and
or r.univar for basic descriptive stats. If all looks fine (min value,
max value, mean and standard deviation show clearly there is useful
information in the queried image(s)), then it's only a color-table
issue. You should only need to apply the correct color table via
r.colors.

   Forunately, Moritz was right, thanks a lot! Now it looks normal.

Sure. But consider what is your aim and if you are ok giving away
some portion of the radiometric details.

   And you wrote also that "the correct way is to not use the Digital Numbers
   directly. Rather, convert to reflectance, then do whatever you have to."
   So application of i.landsat.toar is correct, right?

Yes.

   Sorry for all stupid/basic questions but this is my first serious time
   with grass and landsat image processing.

If you do request (from me, at least) to accept apologies, then accept
my apologies for my infinite silly questions I have asked in this list
:slight_smile:

Nikos

ps- Please try to keep some flow-of-text from top-to-bottom when
replying in the list.

Hello! It’s me again :slight_smile:

If you go all the way from DN to Top-of-Atmosphere reflectances (via

i.landsat.toar), then to Top-of-Canopy (via i.atcorr), you’ll have

floating point values ranging in [0, 1.0]. If you recode this back to

8-bit,

you should consider whether there is an important loss of the

radiometric resolution.

So, recoding to another range is different than converting to integer

numbers and trying to preserve the range.

The thing that worries me is that I don’t know how to check if the loss of the radiometric resolution is important :confused: What should I compare?

I was trying to convert to integer through r.mapcalc (the only other way I’ve found) with the function int(x) and x was my map (I hope it was ok to put map instead of x) but the result was a map with “one shade” of grey, so maybe it was wrong to put a map instead of x…

According to i.atcorr there is an option “output raster map as integer” (i), so if my input file will be _toar2@konfa (float) and I check that option I will have a map with integer values, right?

However, the most confusing thing for me with i.atcorr is “the aerosol optical depth”. I don’t have “meteorological parameter visibility”. I have images from 1984 and 2007. I’ve found files for Global Aerosol Climatology (1981-2006) on this website http://gacp.giss.nasa.gov/data_sets/ I’ve found the proper row and column in the asci format, but they don’t have data for 2007. I was trying to find it on different pages http://modis-atmos.gsfc.nasa.gov/MOD04_L2/acquiring.html but areas of my interest are black - so no data. On your page github with i.landsat.toar and i.atcorr you wrote "the value for aerosols optical depth (AOD), is set to 0.111 for winter and 0.222 for summer aquisitions" Are these default values? For DOS methods?

Can I use i.landsat.toar with DOS3 instead of i.atcorr? The others where using this http://gis.stackexchange.com/questions/126742/which-dos-method-use-to-convert-at-sensor-radiance-to-at-surface-values-in-grass And also on your graph (the one on github) Nikos, this DOS method leads to “corrected” reflectance. So it means that I can omit i.atccor, right? I’m thinking about preprocessing of my images like this: i.landsat.toar + DOS3, then r.recode (I don’t know the other way to change float to integer), then i.histo.match And after that classification (I want to compare land change cover, it is Syrian Jezirah, so +/-95% is agriculture).

Greetings

Joanna

On 05/05/15 17:18, joanna mardas wrote:

Hello! It's me again :slight_smile:

> If you go all the way from DN to Top-of-Atmosphere reflectances (via

> i.landsat.toar), then to Top-of-Canopy (via i.atcorr), you'll have

> floating point values ranging in [0, 1.0]. If you recode this back to

> 8-bit,

> you should consider whether there is an important loss of the

> radiometric resolution.

>

> So, recoding to another range is different than converting to integer

> numbers and trying to preserve the range.

The thing that worries me is that I don't know how to check if the loss
of the radiometric resolution is important :confused: What should I compare?

I'm not sure that for Landsat 5 the loss is so important, but you can visually compare an image recoded to 0-255 with the one coming out of i.landsat.toar...

I was trying to convert to integer through *r.mapcalc* (the only other
way I've found) with the function *int(x)* and x was my map (I hope it
was ok to put map instead of x) but the result was a map with "one
shade" of grey, so maybe it was wrong to put a map instead of x...

int() rounds, so if you maps varies between 0 and 1, rounding will make you lose almost all information. You have to multiply your values by a multiple of ten corresponding to the degree of precision you want before converting it to int. I.e. for a value of 0.12345678, if you want to keep the first five numbers after the decimal point, you have to do something like this:

newMap = int(oldMap * 100000)

According to*i.atcorr* there is an option "output raster map as integer"
(i), so if my input file will be _toar2@konfa (float) and I check that
option I will have a map with integer values, right?

Yes.

However, the most confusing thing for me with i.atcorr is /"the aerosol
optical depth"/. I don't have "meteorological parameter visibility". I
have images from 1984 and 2007. I've found files for Global Aerosol
Climatology (1981-2006) on this
website http://gacp.giss.nasa.gov/data_sets/ I've found the proper row
and column in the asci format, but they don't have data for 2007. I was
trying to find it on different
pages http://modis-atmos.gsfc.nasa.gov/MOD04_L2/acquiring.html but areas
of my interest are black - so no data. On your page github with
i.landsat.toar and i.atcorr you wrote"/the value for aerosols optical
depth (AOD), is set to 0.111 for winter and 0.222 for summer
aquisitions" /Are these default values? For DOS methods?

Yes, these are default values which might not be applicable to your case.

AFAIU, you could use the vibility in km instead of the aerosol optical depth. i.atcorr will then calculate optical depth based on a standard model. To get visibility find a meteorological station close to your image and get the info for the relevant date from that station.

A site I find convenient for that is http://www.ogimet.com/.

Can I use*i.landsat.toar with DOS3 * instead of i.atcorr? The others
where using
this http://gis.stackexchange.com/questions/126742/which-dos-method-use-to-convert-at-sensor-radiance-to-at-surface-values-in-grass
And also on your graph (the one on github) Nikos, this DOS method leads
to "corrected" reflectance. So it means that I can omit i.atccor, right?

i.atcorr provides a much more sophisticated algorithm for atmospheric corrections. This does not necessarily mean that is "more correct", but at least that it tries to take into account more information. DOS is a very simple algorithm of correction with the advantage that you don't need much info to run it, but it is generally even more of an approximation than i.atcorr.

I'm thinking about preprocessing of my images like this:*i.landsat.toar
+ DOS3*, then*r.recode* (I don't know the other way to change float to
integer), then*i.histo.match* And after that classification

I'm not sure that histogram matching is really necessary, or even advisable, before classification, but why not try the path you propose and then check the validity of the results. If they are not good enough for your purpose then you can try to improve by using more sophisticated correction.

Moritz

joanna mardas:

> I was trying to convert to integer through *r.mapcalc* (the only other
> way I've found) with the function *int(x)* and x was my map (I hope it
> was ok to put map instead of x) but the result was a map with "one
> shade" of grey, so maybe it was wrong to put a map instead of x...

Moritz:

int() rounds, so if you maps varies between 0 and 1, rounding will make
you lose almost all information. You have to multiply your values by a
multiple of ten corresponding to the degree of precision you want before
converting it to int. I.e. for a value of 0.12345678, if you want to
keep the first five numbers after the decimal point, you have to do
something like this:

newMap = int(oldMap * 100000)

I wrote on "Mon, 4 May 2015 22:14:41 +0300" how I approached this in the
past. will try to support, if I can, for the rest of questions as well.

Thanks, Nikos

[snip]

Nikos:

> > If you go all the way from DN to Top-of-Atmosphere reflectances (via
> > i.landsat.toar), then to Top-of-Canopy (via i.atcorr), you'll have
> > floating point values ranging in [0, 1.0]. If you recode this back to
> > 8-bit, you should consider whether there is an important loss of the
> > radiometric resolution.
> > So, recoding to another range is different than converting to integer
> > numbers and trying to preserve the range.

joanna:

> The thing that worries me is that I don't know how to check if the loss
> of the radiometric resolution is important :confused: What should I compare?

Moritz:

I'm not sure that for Landsat 5 the loss is so important, but you can
visually compare an image recoded to 0-255 with the one coming out of
i.landsat.toar...

Nor am I sure about it. Landsat5 is 8-bit. But one should definitively consider it, and mention the
decisions taken while documenting the process.

[some text removed]

joanna:

> According to*i.atcorr* there is an option "output raster map as integer"
> (i), so if my input file will be _toar2@konfa (float) and I check that
> option I will have a map with integer values, right?

Moritz:

Yes.

joanna:

> However, the most confusing thing for me with i.atcorr is /"the aerosol
> optical depth"/. I don't have "meteorological parameter visibility". I
> have images from 1984 and 2007. I've found files for Global Aerosol
> Climatology (1981-2006) on this
> website http://gacp.giss.nasa.gov/data_sets/ I've found the proper row
> and column in the asci format, but they don't have data for 2007. I was
> trying to find it on different
> pages http://modis-atmos.gsfc.nasa.gov/MOD04_L2/acquiring.html but areas
> of my interest are black - so no data. On your page github with
> i.landsat.toar and i.atcorr you wrote"/the value for aerosols optical
> depth (AOD), is set to 0.111 for winter and 0.222 for summer
> aquisitions" /Are these default values? For DOS methods?

Moritz:

Yes, these are default values which might not be applicable to your case.

There is a paper, also, suggested by Moritz quite some time ago:
<http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-34-15-2765&gt;\. In
it, there is "Table 4. Rayleigh Optical Depths at 0-km Altitude for Six
Different Atmosphere Models". Perhaps useful.

AFAIU, you could use the vibility in km instead of the aerosol optical
depth. i.atcorr will then calculate optical depth based on a standard
model. To get visibility find a meteorological station close to your
image and get the info for the relevant date from that station.

Yes, that should do as well.

A site I find convenient for that is http://www.ogimet.com/.

Thanks for the tip :slight_smile:

> Can I use*i.landsat.toar with DOS3 * instead of i.atcorr? The others
> where using
> this http://gis.stackexchange.com/questions/126742/which-dos-method-use-to-convert-at-sensor-radiance-to-at-surface-values-in-grass
> And also on your graph (the one on github) Nikos, this DOS method leads
> to "corrected" reflectance. So it means that I can omit i.atccor, right?

i.atcorr provides a much more sophisticated algorithm for atmospheric
corrections. This does not necessarily mean that is "more correct", but
at least that it tries to take into account more information. DOS is a
very simple algorithm of correction with the advantage that you don't
need much info to run it, but it is generally even more of an
approximation than i.atcorr.

> I'm thinking about preprocessing of my images like this:*i.landsat.toar
> + DOS3*, then*r.recode* (I don't know the other way to change float to
> integer), then*i.histo.match* And after that classification

joanna, once again, the easy "other way" is posted in my first reply, I
think. You just need to multiply with 1000, perform the histogram
matching, then divide by 1000.0 to get back to floats.

I'm not sure that histogram matching is really necessary, or even
advisable, before classification, but why not try the path you propose
and then check the validity of the results. If they are not good enough
for your purpose then you can try to improve by using more sophisticated
correction.

As we all know, if one tries to compare scenes over the same area, aqcuired at different
times, it's necessary to relatively normalise'em (different dates,
different solar geometries, variations in the spectral response of the
same surfaces). The same, I think, is valid if one combines multiple
scenes acquired at the same date but cover adjacent areas.
A relative normalisation can help, in such cases, a lot to make classification
results comparable. For the latter, perhaps it is not necessary in flat
areas!?

Of course, if the approach is going to be one, independent,
classification per scene, and then try to compare the outcomes, things
are very different. It might work well without undergoing relative
normalisation actions.

Histogram matching could be used as a mean for relative atmospheric correction.

Also, there is an effort to do something more sophisticated in this direction
by Tomas Brunclik:
<http://www.researchgate.net/publication/275020325_i.grid.correl.atcor_version_0.91b&gt;\.

I haven't checked what's the latest status of it, nor had I any contact
with the author recently (we did discuss something in the past).

I am very interested in his work as I have
performed similar computations in the past using messy scripts in GRASS
combined with some linear regressions in R. Maybe his tool is more
mature now?

Nikos

Thanks Nikos!

> I'm not sure that for Landsat 5 the loss is so important, but you can
> visually compare an image recoded to 0-255 with the one coming out of
> i.landsat.toar...

Nor am I sure about it. Landsat5 is 8-bit. But one should definitively
consider it, and mention the
decisions taken while documenting the process.

I did both, r.recode to 0-255 and r.mapcalc int(oldmap*100000). Both output images and also histograms look ok I think. The histograms have the same shapes and visually images are the same, maybe just slight differences but it's impossible to notice.

There is a paper, also, suggested by Moritz quite some time ago:
<http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-34-15-2765&gt;\. In
it, there is "Table 4. Rayleigh Optical Depths at 0-km Altitude for Six
Different Atmosphere Models". Perhaps useful.

I'll definitelly look at this.

joanna, once again, the easy "other way" is posted in my first reply, I
think. You just need to multiply with 1000, perform the histogram
matching, then divide by 1000.0 to get back to floats.

I'll do DOS3, cuz I have to read more about aerosol depth etc. So for now, I think that DOS will be better. And I'll do the histogram matching. I've tried to do this once (after r.recode) but I think (well I'm sure) that I did this in a wrong way cuz I've matched all bands from both 1984 and 2007. It looked really bad haha :smiley: I guess I should match band to band, for example landsat07B2 to landsat84B2, landsat07B3 to landsat84B3 etc.

As we all know, if one tries to compare scenes over the same area,
aqcuired at different
times, it's necessary to relatively normalise'em (different dates,
different solar geometries, variations in the spectral response of the
same surfaces). The same, I think, is valid if one combines multiple
scenes acquired at the same date but cover adjacent areas.
A relative normalisation can help, in such cases, a lot to make
classification
results comparable. For the latter, perhaps it is not necessary in flat
areas!?

I agree, that it should be done. Well my area is flat almost like a table, there are just "tells" (artificial settlement hills, ancient).

Of course, if the approach is going to be one, independent,
classification per scene, and then try to compare the outcomes, things
are very different. It might work well without undergoing relative
normalisation actions.

That is what I was going to do, an independent classification per scene. I have only two so it's not a big problem.

Histogram matching could be used as a mean for relative atmospheric
correction.

That is interesting.

Also, there is an effort to do something more sophisticated in this
direction
by Tomas Brunclik:
<http://www.researchgate.net/publication/275020325_i.grid.correl.atcor_version_0.91b&gt;\.

I haven't checked what's the latest status of it, nor had I any contact
with the author recently (we did discuss something in the past).

I am very interested in his work as I have
performed similar computations in the past using messy scripts in GRASS
combined with some linear regressions in R. Maybe his tool is more
mature now?

Next time maybe, now I'm too green for this :slight_smile: R... I've heard about this and that's the end of my knowledge about R (sorry, I'm just an archaeologist).
Joanna :slight_smile:

On 07/05/15 20:06, Nikos Alexandris wrote:

Moritz:

AFAIU, you could use the vibility in km instead of the aerosol optical
depth. i.atcorr will then calculate optical depth based on a standard
model. To get visibility find a meteorological station close to your
image and get the info for the relevant date from that station.

Yes, that should do as well.

A site I find convenient for that is http://www.ogimet.com/.

Thanks for the tip :slight_smile:

It might actually be interesting to explore existing web APIs for meteorological data. ogimet provides one[1], but I think there are others (e.g. Meteomanz.com [2]). Metaf2xml also seems to provide access to the different APIs [3].

This might allow to automagically get the visibility data from the nearest station for the date and time closest to image data and time.

Just an idea for someone to work on in the future...

:wink:

Moritz

[1] http://www.ogimet.com/getsynop_help.phtml.en
[2] http://www.meteomanz.com/index?l=1
[3] http://metaf2xml.sourceforge.net/

Hey!

joanna, once again, the easy "other way" is posted in my first reply, I
think. You just need to multiply with 1000, perform the histogram
matching, then divide by 1000.0 to get back to floats.

So I did that: int(oldmap*1000), then i.histo.match (band to band). I did also r.recode (oldlow:oldhigh:0:255), and then i.histo.match (band to band). Histograms in the middle (results of int(oldmap*1000) and i.hist.match) look strange and the images have "holes"-places with no data, totally transparent. After r.recode there are no places with no data or there are 1 or 2 very small (max 5 pixels). I don't know why it is like that. So I will either do the classification after r.recode and
i.histo.match or without r.recode and i.histo.match.
BTW, does classification work with both float and integer?
Joanna

(attachments)

histograms.jpg

* joanna mardas <joanna.mardas@wp.pl> [2015-05-10 02:28:02 +0200]:

Hey!
> joanna, once again, the easy "other way" is posted in my first reply, I
> think. You just need to multiply with 1000, perform the histogram
> matching, then divide by 1000.0 to get back to floats.

So I did that: int(oldmap*1000), then i.histo.match (band to band). I did also r.recode (oldlow:oldhigh:0:255), and then i.histo.match (band to band). Histograms in the middle (results of int(oldmap*1000) and i.hist.match) look strange and the images have "holes"-places with no data, totally transparent. After r.recode there are no places with no data or there are 1 or 2 very small (max 5 pixels). I don't know why it is like that. So I will either do the classification after r.recode and
i.histo.match or without r.recode and i.histo.match.

All right :-). Try what I posted on [2015-05-04 22:14:41]:

# if Reflectance > 1, set to 1000, else round and multiply by 1000
r.mapcalc --o "TopoCorr_B_Trimmed_DOS1_ToCR_${Band_No}_${SCENE} = if(
TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} > 1, 1000, round( 1000
* TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} ) )"

# histogram matching
...

# convert histo-matched images back to floating point: rescale to 0-1.0
r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"

This is:

- don't use the int() function
- just mulitply by 1000 and round, eg:
r.mapcalc "output = if( inputmap > 1, 1000, round( 1000 * inputmap))"

- do the histogram matching
- afterwards, divide by '1000.0' # the dot is important!

BTW, does classification work with both float and integer?
Joanna

I think so. Check the manuals.

Best of success, Nikos

Hey!

All right :-). Try what I posted on [2015-05-04 22:14:41]:

# if Reflectance > 1, set to 1000, else round and multiply by 1000
r.mapcalc --o "TopoCorr_B_Trimmed_DOS1_ToCR_${Band_No}_${SCENE} = if(
TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} > 1, 1000, round( 1000
* TopoCorr.B.Trimmed.DOS1.ToCR.${Band_No}@${SCENE} ) )"

# histogram matching
...

# convert histo-matched images back to floating point: rescale to 0-1.0
r.mapcalc --o "${HistoMatchedMap} = ${HistoMatchedMap} / 1000.0"

This is:

- don't use the int() function
- just mulitply by 1000 and round, eg:
r.mapcalc "output = if( inputmap > 1, 1000, round( 1000 * inputmap))"

- do the histogram matching
- afterwards, divide by '1000.0' # the dot is important!

So I did all of that and there were still holes on images and spikes in the graph... so I did some "experiments" on bands 5 (the biggest holes). I've matched original bands, in result no holes in images, but spikes in the graph. I did again i.landsat.toar + this time DOS1 (cuz I thought that maybe it's the fault of DOS3), and then r.mapcalc and i.histo.match as you've said. After that, there were holes and spikes. I've checked raster statistics (for images after toar+dos and r.mapcalc) and its
max and min values. The image's max value for landsat84 was 248 and for landsat07 346 (after toar+dos and r.mapcalc). The default value for parameter max in i.histo.match is 255, so I've changed it to 346 to see what's gonna happen. And the output images have no holes, spikes in the graphs are still present. But since there are spikes even in graphs of original landsat images, maybe the spikes are normal situation. Maybe the holes and spikes occure when the difference between max values of
images is too big (in case of my bands 5 almost 100), after i.histo.match of r.recoded (to 0:255) images (band 4) there were no holes and no spikes - no or small differences between max values.
Greetings
Joanna