[GRASS-dev] New georectifying module in TclTk

[cmbarton]

I¹ve nearly finished a new georectifying module in TclTk and want to get it
into the CVS where it can be tested and I can get some help on RMS error
calculation.

Here is a brief rundown on how it works.

Start out in the location/mapset where you want the georectified map to END
UP (i.e., the TARGET location). Open any map display (or displays) that you
want to use for getting geographic coordinates

Select georectify from the file menu

In the georectify startup, you must:
    1) select the mapset of the file to be georectified (the location and
gisdbase are automatically noted)
    2) run i.group to create a group to georectify. Target is automatically
set to the current location and mapset.
    3) select the group to georectify
    4) select the raster map to display for interactive georectification
(normally a raster in the group you want to georectify)
    5) start the georectification

You will have a map display with the raster to georectify and a GCP (ground
control point) manager window.
    1) Click in an xy entry box in the GCP manager
    2) Click on the xy map you want to georectify to put xy coordinates in
the entry box. The focus will shift to the corresponding geographic entry
box in the GCP manager.
    3) Click on a point in a georectified map to get geographic coordinates.
    4) Repeat until you have enough xy and geographic coordinate pairs.

    NOTE: you can edit/enter any point by typing. You can also edit
existing points by clicking on a new spot in the display.

Click the button to save the points to a POINTS file

Click on the georectify button to select the type of georectification you
want (Is 3rd order georectification still broken?)

That's it.

RMS error calculation does not work yet. There may be some issues still
remaining with zooming in the xy map display.

***************

Anyone who is familiar with how i.rectify works, I need to ask for help on
RMS error calculation. Somehow, I need to actually rectify the GCP's so I
will get meainingful RMS error results. I can't find where to do this in the
C code (semi-incomprehensible to me).

Try it out and let me know what you think.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics and Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

[Glynn_Clements1]

Michael Barton wrote:

Anyone who is familiar with how i.rectify works, I need to ask for help on
RMS error calculation. Somehow, I need to actually rectify the GCP's so I
will get meainingful RMS error results. I can't find where to do this in the
C code (semi-incomprehensible to me).

AFAICT, the relevant function is compute_transformation() in
imagery/i.points/analyze.c.

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

[cmbarton]

I will admit from the beginning that since I don't really know C, I may well
be missing something. But I think this calls some transformations in
i.rectify somehow. As I understand it (also in the i.points and i.rectify
docs), you first need to rectify the points using one of the rectification
transformations. Then you can go ahead with the compute_transformation()
procedure.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Glynn Clements <glynn@gclements.plus.com>
Date: Tue, 6 Jun 2006 02:17:29 +0100
To: Michael Barton <michael.barton@asu.edu>
Cc: GRASS developers list <grass-dev@grass.itc.it>, Helena Mitasova
<hmitaso@unity.ncsu.edu>
Subject: Re: [GRASS-dev] New georectifying module in TclTk

Michael Barton wrote:

Anyone who is familiar with how i.rectify works, I need to ask for help on
RMS error calculation. Somehow, I need to actually rectify the GCP's so I
will get meainingful RMS error results. I can't find where to do this in the
C code (semi-incomprehensible to me).

AFAICT, the relevant function is compute_transformation() in
imagery/i.points/analyze.c.

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

[Glynn_Clements1]

Michael Barton wrote:

I will admit from the beginning that since I don't really know C, I may well
be missing something. But I think this calls some transformations in
i.rectify somehow. As I understand it (also in the i.points and i.rectify
docs), you first need to rectify the points using one of the rectification
transformations. Then you can go ahead with the compute_transformation()
procedure.

compute_transformation() starts by calling Compute_equation() (in
imagery/i.points/equ.c), which calls I_compute_georef_equations().
That computes the best-fit affine transformation (forward and inverse)
from the control points.

It then goes on to call I_georef() for each pair of control points,
measuring the difference between the projected source and the
destination, and vice-versa.

I_compute_georef_equations() and I_georef() are both defined in
lib/imagery/georef.c.

On entry, the only values which need to have been set are the control
points.

Upon completion:

+ group.equation_stat contains 1 if a transformation could be
computed.

+ group.{E,N}{12,21} contain the transformation coefficients (see
I_georef() for their interpretation).

+ xres, yres and gnd contain the x, y and diagonal (Euclidian)
differences for each pair.

+ xmax, ymax and gmax contain the indices of the control points with
the greatest x, y and diagonal (Euclidian) differences.

+ rms is set to the RMS diagonal error

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

[cmbarton]

Shows how little I understand C. This will help to try and follow out the
various transformations. Mainly I need the affine transformation.

A question here. Is a single transformation used to generate the projected
points for RMS error calcuations? Or are different transformations used
depending on whether 1st, 2nd, or 3rd order is selected?

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Glynn Clements <glynn@gclements.plus.com>
Date: Tue, 6 Jun 2006 05:52:54 +0100
To: Michael Barton <michael.barton@asu.edu>
Cc: Helena Mitasova <hmitaso@unity.ncsu.edu>, GRASS developers list
<grass-dev@grass.itc.it>
Subject: Re: [GRASS-dev] New georectifying module in TclTk

Michael Barton wrote:

I will admit from the beginning that since I don't really know C, I may well
be missing something. But I think this calls some transformations in
i.rectify somehow. As I understand it (also in the i.points and i.rectify
docs), you first need to rectify the points using one of the rectification
transformations. Then you can go ahead with the compute_transformation()
procedure.

compute_transformation() starts by calling Compute_equation() (in
imagery/i.points/equ.c), which calls I_compute_georef_equations().
That computes the best-fit affine transformation (forward and inverse)
from the control points.

It then goes on to call I_georef() for each pair of control points,
measuring the difference between the projected source and the
destination, and vice-versa.

I_compute_georef_equations() and I_georef() are both defined in
lib/imagery/georef.c.

On entry, the only values which need to have been set are the control
points.

Upon completion:

+ group.equation_stat contains 1 if a transformation could be
computed.

+ group.{E,N}{12,21} contain the transformation coefficients (see
I_georef() for their interpretation).

+ xres, yres and gnd contain the x, y and diagonal (Euclidian)
differences for each pair.

+ xmax, ymax and gmax contain the indices of the control points with
the greatest x, y and diagonal (Euclidian) differences.

+ rms is set to the RMS diagonal error

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

[cmbarton]

Where is this procedure ( I_compute_georef_equations())? It isn't in equ.c.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Glynn Clements <glynn@gclements.plus.com>
Date: Tue, 6 Jun 2006 05:52:54 +0100
To: Michael Barton <michael.barton@asu.edu>
Cc: Helena Mitasova <hmitaso@unity.ncsu.edu>, GRASS developers list
<grass-dev@grass.itc.it>
Subject: Re: [GRASS-dev] New georectifying module in TclTk

which calls I_compute_georef_equations().
That computes the best-fit affine transformation (forward and inverse)
from the control points.

[Maciej_Sieczka3]

On Mon, 05 Jun 2006 17:28:03 -0700
Michael Barton <michael.barton@asu.edu> wrote:

<snip>

Congrats Michael! I'm not able to do any testing but will do it in
future and let you know how it works.

Click on the georectify button to select the type of georectification
you want (Is 3rd order georectification still broken?)

Broken?! Can you elaborate on it? The last time I used i.rectify few
mnonths ago I only noticed that 2nd order output is somewhat strange
(map borders became convex, in general rectification result worse
than 1st and 3rd order), but 1st and 3rd produced a reasonable output...

Maciek

--------------------
W polskim Internecie s? setki milion?w stron. My przekazujemy Tobie tylko najlepsze z nich!
http://katalog.panoramainternetu.pl/

[Glynn_Clements1]

Michael Barton wrote:

Shows how little I understand C. This will help to try and follow out the
various transformations. Mainly I need the affine transformation.

A question here. Is a single transformation used to generate the projected
points for RMS error calcuations? Or are different transformations used
depending on whether 1st, 2nd, or 3rd order is selected?

i.points only understands affine (first-order) transformations, so
that's what the error estimates assume.

The 2nd- and 3rd-order transformations are specific to i.rectify, in
imagery/i.rectify/crs.c.

The problem with higher-order transformations is that you can't
readily invert them, so you can only calculate the error in one
direction.

With regard to your subsequent message:

> I_compute_georef_equations() and I_georef() are both defined in
> lib/imagery/georef.c.

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

[H_B]

Where is this procedure ( I_compute_georef_equations())? It isn't in
equ.c.

I_*() are in lib/imagery/
G_*() are in lib/gis/
R_*() are in lib/raster/
etc

$ cd lib/
$ grep -rI "I_compute_georef_equations" *
imagery/georef.c:int I_compute_georef_equations(

so it is in lib/imagery/georef.c

Hamish

[H_B]

Michael:

> (Is 3rd order georectification still broken?)

Maciek:

Broken?! Can you elaborate on it? The last time I used i.rectify few
mnonths ago I only noticed that 2nd order output is somewhat strange
(map borders became convex, in general rectification result worse
than 1st and 3rd order), but 1st and 3rd produced a reasonable
output...

see
https://intevation.de/rt/webrt?serial_num=3166
related ?
https://intevation.de/rt/webrt?serial_num=3052

also
https://intevation.de/rt/webrt?serial_num=3296

I think GDAL order=3 will have the same problems, as both started from
the same code.

I didn't know order=2 had problems..?

Michael:

> There has been some discussion about dropping i.rectify for
> gdalwarp. I¹ve used i.rectify here, but it is probably possible to
> rewrite this to use gdalwarp in the future.

I think we should keep i.rectify but use GDAL's warp API for the heavy
lifting. (as opposed to the "gdalwarp" program)

see
https://intevation.de/rt/webrt?serial_num=2952
http://www.gdal.org/warptut.html

Hamish

[cmbarton]

Thanks. That is helpful. It will be certainly easier if I only am doing rms
calculations using one equation. Still, I wonder can minimizing errors using
an affine transformation produce poorer results for polynomial
transformations?

Also, in responding to my own comment yesterday... I looked at gdalwarp last
night and I'm not sure it could replace i.rectify as is. If there is someway
to specify a set of GCP's that are not embedded in the raster file format, I
didn't see it in the docs. Also, I didn't see any way of using it to do
error calculations (though this might be moot). GDAL will also need to to
GRASS to GRASS file formats. I know this was a problem awhile back, but
perhaps it is not one now.

More interestingly, it might be more possible to wrap v.transform into the
georectifier so that it will do both raster and vector files, if this is of
any interest to people.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics & Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Glynn Clements <glynn@gclements.plus.com>
Date: Tue, 6 Jun 2006 08:46:29 +0100
To: Michael Barton <Michael.Barton@asu.edu>
Cc: Helena Mitasova <hmitaso@unity.ncsu.edu>, GRASS developers list
<grass-dev@grass.itc.it>
Subject: Re: [GRASS-dev] New georectifying module in TclTk

Michael Barton wrote:

Shows how little I understand C. This will help to try and follow out the
various transformations. Mainly I need the affine transformation.

A question here. Is a single transformation used to generate the projected
points for RMS error calcuations? Or are different transformations used
depending on whether 1st, 2nd, or 3rd order is selected?

i.points only understands affine (first-order) transformations, so
that's what the error estimates assume.

The 2nd- and 3rd-order transformations are specific to i.rectify, in
imagery/i.rectify/crs.c.

The problem with higher-order transformations is that you can't
readily invert them, so you can only calculate the error in one
direction.

With regard to your subsequent message:

I_compute_georef_equations() and I_georef() are both defined in
lib/imagery/georef.c.

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

[Maciej_Sieczka3]

On Tue, 06 Jun 2006 08:28:14 -0700
Michael Barton <michael.barton@asu.edu> wrote:

Thanks. That is helpful. It will be certainly easier if I only am
doing rms calculations using one equation. Still, I wonder can
minimizing errors using an affine transformation produce poorer
results for polynomial transformations?

Also, in responding to my own comment yesterday... I looked at
gdalwarp last night and I'm not sure it could replace i.rectify as
is. If there is someway to specify a set of GCP's that are not
embedded in the raster file format, I didn't see it in the docs.

1. 'gdal_translate -gcp' to add GCPs to raster
2. 'gdalwarp' to warp using those GCPs

Also, I didn't see any way of using it to do error calculations
(though this might be moot).

Not in the verbose gdalwarp output AFAIK.

GDAL will also need to to GRASS to GRASS
file formats. I know this was a problem awhile back, but perhaps it
is not one now.

AFAIK Grass raster and vectors are supported only for reading. I've
been using it for e.g. QGIS and it works very good.

More interestingly, it might be more possible to wrap v.transform
into the georectifier so that it will do both raster and vector
files, if this is of any interest to people.

Maybe a good idea. If this also enables us to overlay raster and vector
in the reference window, to select GCPs for raster from vectors and
vice versa, it could be really usefull.

Maciek

--------------------
W polskim Internecie s? setki milion?w stron. My przekazujemy Tobie tylko najlepsze z nich!
http://katalog.panoramainternetu.pl/

[cmbarton]

Hamish,

It sounds like most of the issues can be avoided by using i.rectify -c. What
are the disadvantages of using the -c switch? It is not very clear (to me)
in the docs.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics and Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Hamish <hamish_nospam@yahoo.com>
Date: Tue, 6 Jun 2006 22:35:32 +1200
To: Maciek Sieczka <werchowyna@epf.pl>
Cc: <michael.barton@asu.edu>, <grass-dev@grass.itc.it>
Subject: Re: [GRASS-dev] New georectifying module in TclTk

Michael:

(Is 3rd order georectification still broken?)

Maciek:

Broken?! Can you elaborate on it? The last time I used i.rectify few
mnonths ago I only noticed that 2nd order output is somewhat strange
(map borders became convex, in general rectification result worse
than 1st and 3rd order), but 1st and 3rd produced a reasonable
output...

see
https://intevation.de/rt/webrt?serial_num=3166
related ?
https://intevation.de/rt/webrt?serial_num=3052

also
https://intevation.de/rt/webrt?serial_num=3296

I think GDAL order=3 will have the same problems, as both started from
the same code.

I didn't know order=2 had problems..?

Michael:

There has been some discussion about dropping i.rectify for
gdalwarp. I¹ve used i.rectify here, but it is probably possible to
rewrite this to use gdalwarp in the future.

I think we should keep i.rectify but use GDAL's warp API for the heavy
lifting. (as opposed to the "gdalwarp" program)

see
https://intevation.de/rt/webrt?serial_num=2952
http://www.gdal.org/warptut.html

Hamish

[cmbarton]

For RMS error calculation, I just need to use either the forward or reverse,
right? I don't see where I need both directions.

Michael
__________________________________________
Michael Barton, Professor of Anthropology
School of Human Evolution & Social Change
Center for Social Dynamics and Complexity
Arizona State University

phone: 480-965-6213
fax: 480-965-7671
www: http://www.public.asu.edu/~cmbarton

From: Glynn Clements <glynn@gclements.plus.com>
Date: Tue, 6 Jun 2006 05:52:54 +0100
To: Michael Barton <michael.barton@asu.edu>
Cc: Helena Mitasova <hmitaso@unity.ncsu.edu>, GRASS developers list
<grass-dev@grass.itc.it>
Subject: Re: [GRASS-dev] New georectifying module in TclTk

compute_transformation() starts by calling Compute_equation() (in
imagery/i.points/equ.c), which calls I_compute_georef_equations().
That computes the best-fit affine transformation (forward and inverse)
from the control points.

[Glynn_Clements1]

Michael Barton wrote:

Thanks. That is helpful. It will be certainly easier if I only am doing rms
calculations using one equation. Still, I wonder can minimizing errors using
an affine transformation produce poorer results for polynomial
transformations?

If you intend to use a polynomial transformation, the error
measurements for an affine transformation aren't really relevant.

In any case, the error measurement is only really useful to

Note that the fitting and measurement performed by i.points has no
effect upon the transformation chosen by i.rectify. It's just giving
the user an idea of what to expect from i.rectify (assuming that an
affine transformation is used), as well as highlighting points which
the user might have misplaced.

Also, in responding to my own comment yesterday... I looked at gdalwarp last
night and I'm not sure it could replace i.rectify as is. If there is someway
to specify a set of GCP's that are not embedded in the raster file format, I
didn't see it in the docs. Also, I didn't see any way of using it to do
error calculations (though this might be moot). GDAL will also need to to
GRASS to GRASS file formats. I know this was a problem awhile back, but
perhaps it is not one now.

It would be useful if the same front-end could be used for both
i.rectify and gdalwarp. The main issue is that you need a way to
compute the best-fit transformation from the control points and to
project sample points without actually transformating the entire
image.

More interestingly, it might be more possible to wrap v.transform into the
georectifier so that it will do both raster and vector files, if this is of
any interest to people.

v.transform also uses an affine transform. Non-affine transformations
on vectors are awkward as they transform straight line segments to
curve segments. Also, vector transformations are normally done in the
opposite direction to raster transformations.

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

[Glynn_Clements1]

Michael Barton wrote:

For RMS error calculation, I just need to use either the forward or reverse,
right? I don't see where I need both directions.

You probably don't really need both, at least for an affine
transformation. The diagonal error computed in one direction will be
proportional to that computed in the other direction. The same doesn't
apply for the x and y error, though.

For polynomial transformations, the situation is more complex. The
inverse of a polynomial transformation isn't itself a polynomial
transformation. Or, to put it another way, the best-fit reverse
transformation won't be the inverse of the best-fit forward
transformation.

At the very least, you need to make sure to use the correct direction
(rasters are usually transformed by projecting coordinates from the
target location into the image's X/Y coordinate system; at least,
that's how r.proj works).

It might help to compute the error both ways to provide an idea of how
good a fit the transformation is. A transformation which accurately
maps a minimal set of control points won't necessarily accurately map
the entire region.

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

[H_B]

It sounds like most of the issues can be avoided by using i.rectify
-c.

that is my understanding.

What are the disadvantages of using the -c switch?

you have to restart GRASS in the target location and set up the region
by hand. Ok if you knew to project a v.in.region box first and can
calculate your target resolution, & give g.region all the right answers;
but an automatic determination would be much better for most people.

It is not very clear (to me) in the docs.

As it wasn't (isn't) very clear to the folks writing the docs....

You'll have to dig into the archives to read the "region cropping"
history and the debate on the correct way to describe that flag.

AFAIR, the default mode tries not to reproject cells which are not being
used. (off the edge of the page in one projection) That is a very rough
recollection though- it wasn't all to clear to me at the time either.

Hamish

[Maciej_Sieczka3]

On Wed, 7 Jun 2006 21:51:27 +1200
Hamish <hamish_nospam@yahoo.com> wrote:

> It sounds like most of the issues can be avoided by using i.rectify
> -c.

that is my understanding.

> What are the disadvantages of using the -c switch?

you have to restart GRASS in the target location and set up the region
by hand. Ok if you knew to project a v.in.region box first and can
calculate your target resolution

Could you say how to calculate it?

Maciek

--------------------
W polskim Internecie s? setki milion?w stron. My przekazujemy Tobie tylko najlepsze z nich!
http://katalog.panoramainternetu.pl/

[H_B]

Maciek Sieczka wrote:

> > It sounds like most of the issues can be avoided by using
> > i.rectify -c.
>
> that is my understanding.
>
> > What are the disadvantages of using the -c switch?
>
> you have to restart GRASS in the target location and set up the
> region by hand. Ok if you knew to project a v.in.region box first
> and can calculate your target resolution

Could you say how to calculate it?

1) XY to meters:
see "Introduction to GRASS software" by Markus Neteler, 1998 2nd ed.,
part VI, section 1.1.1 (p.27)

2) lat-lon to meters:
1852m per nautical mile
60 nautical miles per degree lat
width of degree lon goes as cos(lat)

so
dist_m_ns = dist_deg * 1852*60
dist_m_ew = dist_deg * 1852*60 * cos(lat)

and

res_ns = dist/rows
res_ew = dist/cols

see also "g.region -e"

Hamish