GRASS GIS: Ticket #929: g.transform: dump coeffs and transform sparse points
https://trac.osgeo.org/grass/ticket/929
<p>
Hi,
</p>
<p>
the attached patch adds a new flag to dump transform coeffs, transform arbitrary input coords to the other coord system, and a adds a flag to tell it to do these things in reverse.
</p>
<p>
it also fixes a little <code>continue</code> bug if you selected format="".
</p>
<p>
(
based on Glynn's tips on the -dev ML.
</p>
<blockquote>
<p>
<a class="ext-link" href="http://thread.gmane.org/gmane.comp.gis.grass.devel/38099"><span class="icon"></span>http://thread.gmane.org/gmane.comp.gis.grass.devel/38099</a>
</p>
</blockquote>
<p>
)
</p>
<p>
I'm not decided what to do about reading from stdin. I suspect I'll have to drop the east_north=x1,y1,x2,y2,xn,yn,... option in favor of an input file containing points to be transformed because the parser insists on two values, and "-,-" if given is read as -0,-0. Well, I guess I could have strcmp() test for "-,-" but it's a bit silly.
I'd rater not add another flag to direct the input.
any ideas?
</p>
<p>
FIXME: if using the reverse flag check the TARGET location's G_projection() type and use it for G_scan_easting() instead of the current location's type.
</p>
<p>
Hamish
</p>
en-usGRASS GIShttps://trac.osgeo.org/grass/chrome/site/grasslogo_vector_small.png
https://trac.osgeo.org/grass/ticket/929
Trac 1.4.3hamishSat, 13 Feb 2010 16:19:07 GMTattachment set
https://trac.osgeo.org/grass/ticket/929
https://trac.osgeo.org/grass/ticket/929
<ul>
<li><strong>attachment</strong>
→ <span class="trac-field-new">gtrans_pts_dmp.diff</span>
</li>
</ul>
<p>
simplify for loops
</p>
TickethamishSat, 13 Feb 2010 16:34:36 GMT
https://trac.osgeo.org/grass/ticket/929#comment:1
https://trac.osgeo.org/grass/ticket/929#comment:1
<p>
(the original idea here was to find the bounding box coords in the target location given the current set of GCPs and region settings)
</p>
TicketglynnSun, 14 Feb 2010 01:57:14 GMT
https://trac.osgeo.org/grass/ticket/929#comment:2
https://trac.osgeo.org/grass/ticket/929#comment:2
<p>
Replying to <a class="closed ticket" href="https://trac.osgeo.org/grass/ticket/929" title="#929: enhancement: g.transform: dump coeffs and transform sparse points (closed: fixed)">hamish</a>:
</p>
<blockquote class="citation">
<p>
I'm not decided what to do about reading from stdin. I suspect I'll have to drop the east_north=x1,y1,x2,y2,xn,yn,... option in favor of an input file containing points to be transformed because the parser insists on two values, and "-,-" if given is read as -0,-0. Well, I guess I could have strcmp() test for "-,-" but it's a bit silly.
I'd rather not add another flag to direct the input.
any ideas?
</p>
</blockquote>
<p>
Don't try and make one option do different things; use a separate option for reading from a file (use file=- for stdin rather than a flag).
</p>
TickethamishSun, 14 Feb 2010 07:08:04 GMTstatus changed; resolution set
https://trac.osgeo.org/grass/ticket/929#comment:3
https://trac.osgeo.org/grass/ticket/929#comment:3
<ul>
<li><strong>status</strong>
<span class="trac-field-old">new</span> → <span class="trac-field-new">closed</span>
</li>
<li><strong>resolution</strong>
→ <span class="trac-field-new">fixed</span>
</li>
</ul>
<p>
cleaned up version committed to 6.5svn and trunk as <a class="changeset" href="https://trac.osgeo.org/grass/changeset/40988" title="add new options to dump coeffs and transform arbitrary points, both ...">r40988</a> and <a class="changeset" href="https://trac.osgeo.org/grass/changeset/40990" title="add new options to dump coeffs and transform arbitrary points, both ...">r40990</a>.
</p>
TickethamishSun, 14 Feb 2010 08:26:06 GMT
https://trac.osgeo.org/grass/ticket/929#comment:4
https://trac.osgeo.org/grass/ticket/929#comment:4
<p>
I have been pointed to this post by the author:
<a class="ext-link" href="http://www.cruisersforum.com/forums/f134/charts-31346-36.html#post398493"><span class="icon"></span>http://www.cruisersforum.com/forums/f134/charts-31346-36.html#post398493</a>
</p>
<p>
with GPL3 code linked at the end of the post. (remove the .doc; stupid webCMS...)
</p>
<p>
quote:
</p>
<hr />
<ul><li>First for all ref points it converts ref lat/lon in ref E/N using the exact projection formulas (Mercator and Transverse-Mercator supported)
</li><li>Second it tries to approximate the ref X/Y to E/N by polynomial approx. It tests 23 polynomials to find the best one (minimum error, well conditioning). It works with at least 2 ref points.
</li><li>Third it uses this couple of transformations (lat/lon->E/N->X/Y) to recalculate 36 new ref points with a disposal that is always well-conditioned so it can calculate the final third order polys that goes directly form lat/lon to X/Y and viceversa"
</li></ul><hr />
<p>
/endquote
</p>
<p>
the tests are:
</p>
<pre class="wiki">/*
PolyEleEval - by Marco Certelli - License: GPL - Version 01beta
This function evaluates the elements of the polynomial that approximates a set
of ref points placed in the plane xy. This is specific for XY<->EN transformations
needed in Chart Georeferencing. One approximating poly is defined as:
N = a[0] + a[1]X + a[2]Y + a[3]X^2 + a[4]XY + a[5]Y^2 + a[6]X^3 + a[7]X^2Y + a[8]XY^2 + a[9]Y^3
E = b[0] + b[1]X + b[2]Y + b[3]X^2 + b[4]XY + b[5]Y^2 + b[6]X^3 + b[7]X^2Y + b[8]XY^2 + b[9]Y^3
This function takes the int nref, representing number of available ref points,
2 vectors of double x[] and y[], representing the ref points position (either
x, y or lat, lon), and returns a vector of int e[10] which i-th value [1 or 0]
states if the relative polynomial coefficients a[i] is well or ill-conditioned.
The x[] and y[] input values shall be between 0.0 and 1.0 (the caller shall
scale data). The function itself returns an int 1 if an error occurred.
*/
#define N_POLY_TEST 23
{
int i,j,n,t,ne,optt;
double *A[2*MAX_REF_POINTS+2], B1[2*MAX_REF_POINTS+2], B2[MAX_REF_POINTS], C1[20], C2[10];
double eqmN,eqmE,sN,sE,rcond,sC,prjerror,ill,opt,minopt;
// n X Y X X Y X X X Y a
// X Y Y X X Y Y =
// X Y Y Y b
int test[N_POLY_TEST][11] =
{3,1,1,0,0,0,0,0,0,0,1, // 0 - P=C + X + Y & (Cxx=Cyy, Cxy=-Cyx)
3,1,1,0,0,0,0,0,0,0,0, // 1 - P=C + X + Y
4,1,1,1,0,0,0,0,0,0,1, // 2 - P=C + X + Y + X^2 & (Cxx=Cyy, Cxy=-Cyx)
4,1,1,0,1,0,0,0,0,0,1, // 3 - P=C + X + Y + XY & (Cxx=Cyy, Cxy=-Cyx)
4,1,1,0,0,1,0,0,0,0,1, // 4 - P=C + X + Y + Y^2 & (Cxx=Cyy, Cxy=-Cyx)
4,1,1,1,0,0,0,0,0,0,0, // 5 - P=C + X + Y + X^2
4,1,1,0,1,0,0,0,0,0,0, // 6 - P=C + X + Y + XY
4,1,1,0,0,1,0,0,0,0,0, // 7 - P=C + X + Y + Y^2
5,1,1,1,1,0,0,0,0,0,0, // 8 - P=C + X + Y + X^2 + XY
5,1,1,0,1,1,0,0,0,0,0, // 9 - P=C + X + Y + XY + Y^2
5,1,1,1,0,1,0,0,0,0,0, // 10 - P=C + X + Y + X^2 + Y^2
6,1,1,1,1,1,0,0,0,0,0, // 11 - P=C + X + Y + X^2 + XY + Y^2
6,1,1,1,1,0,0,1,0,0,0, // 12 - P=C + X + Y + X^2 + XY + X^2Y
6,1,1,0,1,1,0,0,1,0,0, // 13 - P=C + X + Y + XY + Y^2+ Y^2X
7,1,1,1,1,1,1,0,0,0,0, // 14 - P=C + X + Y + X^2 + XY + Y^2 + X^3
7,1,1,1,1,1,0,1,0,0,0, // 15 - P=C + X + Y + X^2 + XY + Y^2 + X^2Y
7,1,1,1,1,1,0,0,1,0,0, // 16 - P=C + X + Y + X^2 + XY + Y^2 + XY^2
7,1,1,1,1,1,0,0,0,1,0, // 17 - P=C + X + Y + X^2 + XY + Y^2 + Y^3
8,1,1,1,1,1,0,1,1,0,0, // 18 - P=C + X + Y + X^2 + XY + Y^2 + X^2Y+ XY^2
8,1,1,1,1,1,1,0,0,1,0, // 19 - P=C + X + Y + X^2 + XY + Y^2 + X^3 + Y^3
9,1,1,1,1,1,1,1,1,0,0, // 20 - P=C + X + Y + X^2 + XY + Y^2 + X^3 + X^2Y+ XY^2
9,1,1,1,1,1,0,1,1,1,0, // 21 - P=C + X + Y + X^2 + XY + Y^2 + X^2Y+ XY^2+ Y^3
10,1,1,1,1,1,1,1,1,1,0};// 22 - P=C + X + Y + X^2 + XY + Y^2 + X^3 + X^2Y+ XY^2 + Y^3
int element[N_POLY_TEST][10];
</pre><p>
I am wondering if this is useful for us but I am reminded that in stats that the simplest model fit that explains the data is usually the best one to use. You can use high nth-order polynomials to exactly fit the data, but what you are really doing is fitting the noise.
</p>
<p>
So for example I have a good & undistorted scan of a printed map which I wish to georectify. I can get a smaller RMS by using order=3, but that's just fitting my measurement error better. It's more appropriate to use order=2 which will average out my measurement error instead of modelling the noise -- which is what I <em>really</em> want.
</p>
<p>
comments?
</p>
<p>
Hamish
</p>
TicketglynnSun, 14 Feb 2010 18:38:24 GMT
https://trac.osgeo.org/grass/ticket/929#comment:5
https://trac.osgeo.org/grass/ticket/929#comment:5
<p>
Replying to <a class="ticket" href="https://trac.osgeo.org/grass/ticket/929#comment:4" title="Comment 4">hamish</a>:
</p>
<blockquote class="citation">
<p>
I am wondering if this is useful for us but I am reminded that in stats that the simplest model fit that explains the data is usually the best one to use. You can use high nth-order polynomials to exactly fit the data, but what you are really doing is fitting the noise.
</p>
</blockquote>
<p>
I would expect that you're typically trying to find a suitable approximation to a complex analytic function (i.e. a typical cartographic or geometric projection).
</p>
<p>
For many projections, a quadratic polynomial will be a poor fit. The main reason for using cubic functions for interpolation is that you can specify both the value and first derivative at each end (4 degrees of freedom).
</p>
<p>
If you consider approximating sin(x) symmetrically about zero: you can get a reasonable approximation with a cubic function, but the optimal quadratic approximation will actually be linear.
</p>
<p>
Note that the most complex case in the above code is the same as g.transform with order=3, i.e. generalised cubic (rather than bicubic).
</p>
TickethamishMon, 15 Feb 2010 08:26:34 GMT
https://trac.osgeo.org/grass/ticket/929#comment:6
https://trac.osgeo.org/grass/ticket/929#comment:6
<p>
Replying to <a class="ticket" href="https://trac.osgeo.org/grass/ticket/929#comment:5" title="Comment 5">glynn</a>:
</p>
<blockquote class="citation">
<p>
I would expect that you're typically trying to find a suitable
approximation to a complex analytic function (i.e. a typical
cartographic or geometric projection).
<br /><br />
For many projections, a quadratic polynomial will be a poor fit.
The main reason for using cubic functions for interpolation is
that you can specify both the value and first derivative at each
end (4 degrees of freedom).
</p>
</blockquote>
<p>
Sure, if I were reprojecting or suspected the scan was distorted order=3 would be appropriate-
</p>
<p>
I guess I should have mentioned that I always try very hard to match the target location's projection to the scanned-map's native projection. So I'm using i.rectify mainly as a georeferencing tool, not a reprojection tool. (so rotation+scale but no warping)
</p>
<p>
e.g. if the quality-scan is a nautical chart in Mercator projection, with a known +lat_ts and grid lines in lat/lon, I'll put GCPs on all the grid line confluences* and enter the coords as degrees, and finally run a little script to reproject the target POINTS to the custom merc projection (see wiki addons).
</p>
<p>
Then run r.proj to pull it from the custom merc location into whatever working location it is needed.
</p>
<p>
Probably it is overkill, but trying to do the right thing..
</p>
<p>
Hamish
</p>
<p>
[*] totally unrelated but interesting: <a class="ext-link" href="http://confluence.org/"><span class="icon"></span>http://confluence.org/</a>
</p>
Ticket