I think the basic idea presented by Dozier sounds feasible but his execution
falls short of the goal. Although the the method of computing elliptic integral of the second kind with complex argument may be OK it has troubles with large arguments when evaluating the Jacobian Zeta function. The second, and I believe the more serious problem, is the Newton-Raphson method employed. He has expanded the basic real function and applied it to a complex variable. I am not sure that this is appropriate and searching through the net and all leads me to the conclusion that we are getting into deep water when dealing with the complex plane. Also, it looks like we are also dealing with multiple roots, especially when longitude exceeds a certain value (suggested to be (pi/2)*(1-k))---a factor not addressed in Dozier's solution. Note that the cusp (poorly displayed in the previously mention gif url) appears to be the beginning of the multiple root solution. I say poorly displayed as the gradation of the equitorial parallel should *smoothly* begin a swing to the north OR (importantly) to the south. Selecting the north or south root becomes a practical problem for a projection program and is a problem with any method dealing with the comprehensive TM projection. Because I have no training and no experience in working with complex variable problems and have failed to find any practical material related to alternate methods to compute elliptic integrals with complex arguments and, more importantly, a Newton-Raphson routine for determining roots of complex functions, I have decide to suspend again any activity on the Dozier method. One has to know when to throw in the towel. ;-) -- Jerry and the low-riders: Daisy Mae and Joshua "Cogito cogito ergo cogito sum" Ambrose Bierce, The Devil's Dictionary _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
In a message dated 6/29/06 08:45:00, [hidden email] writes: (Full text of message at end.) ...is the Newton-Raphson method employed. He has expanded the I haven't looked at Dozier in depth, but in general there is nothing less appropriate about Newton-Raphson when applied to complex variables than there is when applied to real-valued variables. It's quite easy to get yourself into trouble when finding roots even of real-valued functions. While the function in question has problematic regions, that has nothing to do with the fact that it is complex-valued. plane. Also, it looks like we are also dealing with multiple roots, I do not understand why you think the fact of multiple roots supports your notion of the intractability of the solution.
The cusp is poorly observed, not poorly displayed. The illustration is exact to sub-pixel resolution and it matches Lee's illustration precisely despite using a different method of calculation. It seems you do not know a continuously differentiable curve when you see one. You have remarked previously that the projection is "not intuitive". To you, of course. Not necessarily to the reading audience. You didn't like the cusp; you thought it indicated the projection wasn't really conformal. Once the cusp was demonstrated (by means of a 30-year old peer-reviewed journal article) to be an attribute of the projection, you decided you didn't like how it looked in an image I supplied, even though that image precisely matches the one on the peer-reviewed article. It is curious that you believe it appropriate to cast public aspersion based on your own, flawed notions rather than objective facts. If that is how you talk yourself out of a project then I don't suppose there is much anyone can do about it, since reason clearly has nothing to do with it. To those interested in the full-spheroid transverse Mercator, I urge you not to be skeptical of its existence or attainability based on Mr. Evenden's comments. There is nothing controversial about either. The mathematics has been published in peer-reviewed journals, confirmed any number of times by people who understand the mathematics, and expressed by at least three different calculational methods (whatever method Lee used; Dozier; and Wallis) in at least five implementations that I know of. While there is treacherous calculational territory to traverse, that is true of many projections. As always, you must understand the domain and choose numerical techniques appropriate to it. To that end, please feel free to contine the discussion on the "Complex Transverse Mercator" thread. Regards, -- daan Strebe In a message dated 6/29/06 08:45:00, [hidden email] writes: I think the basic idea presented by Dozier sounds feasible but his execution _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
On Thursday 29 June 2006 3:30 pm, [hidden email] wrote:
... > > Note that the cusp (poorly displayed in the previously mention gif url) > > appears to be the beginning of the multiple root solution. I say poorly > > displayed as the gradation of the equitorial parallel should *smoothly* > > begin > > a swing to the north OR (importantly) to the south. > > The cusp is poorly observed, not poorly displayed. The illustration is > exact to sub-pixel resolution and it matches Lee's illustration precisely > despite using a different method of calculation. It seems you do not know a > continuously differentiable curve when you see one. My goodness, such a cranky old man! I know you don't like me but I have tried to be civil. Can we bury the hatchet sometime. I am getting too old for this bullshit. -- Jerry and the low-riders: Daisy Mae and Joshua "Cogito cogito ergo cogito sum" Ambrose Bierce, The Devil's Dictionary _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
In reply to this post by strebe
On Thursday 29 June 2006 3:30 pm, [hidden email] wrote:
> In a message dated 6/29/06 08:45:00, [hidden email] writes: > (Full text of message at end.) > > > ...is the Newton-Raphson method employed. He has expanded the > > basic real function and applied it to a complex variable. I am not sure > > that > > this is appropriate... > > I haven't looked at Dozier in depth, but in general there is nothing less > appropriate about Newton-Raphson when applied to complex variables than > there is when applied to real-valued variables. It's quite easy to get > yourself into trouble when finding roots even of real-valued functions. > While the function in question has problematic regions, that has nothing to > do with the fact that it is complex-valued. Isn't that what I just said??? At this point I am picking on Dozier, not you nor anyone else. > > plane. Also, it looks like we are also dealing with multiple roots, > > especially when longitude exceeds a certain value (suggested to be > > (pi/2)*(1-k))---a factor not addressed in Dozier's solution > > I do not understand why you think the fact of multiple roots supports your > notion of the intractability of the solution. Where did I say the solution did not exist nor the plot was not basically correct---just poorly done. ... material previously commented on > You have remarked previously that the projection is "not intuitive". To > you, of course. Not necessarily to the reading audience. You didn't like Let us stop for a moment on this "intuitive issue." Plain old Mercator goes to infinity in the NS directions for both the sphere and elliptical case. When we flip the cylinder over on its side the projection goes to infinity in the east-west direction (as one might logically expect) for the spherical case. SO, is it not intuitive that that the elliptical case would behave in a similar manner? That is all I am saying. I am perfectly willing to accept that for some weird magic of complex analysis that it defies intuition and is finite in all directions. I have accepted this finiteness on faith for several years, but being a good atheist, I must allways view it with skepticism until it is appropriately demonstrated. > the cusp; you thought it indicated the projection wasn't really conformal. > Once the cusp was demonstrated (by means of a 30-year old peer-reviewed > journal article) to be an attribute of the projection, you decided you I had the impression you got the plot from Wallis. Was he peer reviewed? If so, what was the 30 year old publication. > didn't like how it looked in an image I supplied, even though that image > precisely matches the one on the peer-reviewed article. It is curious that > you believe it appropriate to cast public aspersion based on your own, > flawed notions rather than objective facts. If that is how you talk > yourself out of a project then I don't suppose there is much anyone can do > about it, since reason clearly has nothing to do with it. > > To those interested in the full-spheroid transverse Mercator, I urge you > not to be skeptical of its existence or attainability based on Mr. > Evenden's comments. There is nothing controversial about either. The > mathematics has been published in peer-reviewed journals, confirmed any > number of times by people who understand the mathematics, and expressed by > at least three different calculational methods (whatever method Lee used; > Dozier; and Wallis) in at least five implementations that I know of. While > there is treacherous calculational territory to traverse, that is true of > many projections. As always, you must understand the domain and choose > numerical techniques appropriate to it. Stating it another way, my surrender was to the Dozier article that while the method sounds interesting the execution falls short. The code had several errors and some methodology is in question. I have spent an additional week rooting out a complex Bulirsch routine from the old Bell lab site and spent several days relating Jacobi's form to Legendre in order to test against Abramowitz tables and GSL software for at least real arguments---they now agree. But Dozier's code does not agree to better than the third or forth place (I am not saying Dozier's code is wrong yet). And I still have no handle on a replacement for Newton-Raphson. There's the old saying: "I'm from Missouri. Show me." Has someone have a reasonable reference list on this subject (other than Dozier and Wallis). I am aware of Lee but I have not seen some of the other names mentioned without details. > To that end, please feel free to contine the discussion on the "Complex > Transverse Mercator" thread. > > Regards, > -- daan Strebe Lastly, I am not a grad student dependent upon the graces of the almighty professor and thus must speak in muted and humble tones in his presence. Rather than complaining about my skepticism, I'd appreciate a good pointer or hints on complex root finding. :-) -- Jerry and the low-riders: Daisy Mae and Joshua "Cogito cogito ergo cogito sum" Ambrose Bierce, The Devil's Dictionary _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
> From: "Gerald I. Evenden"
> Stating it another way, my surrender was to the Dozier article that while the > method sounds interesting the execution falls short. The code had several > errors and some methodology is in question. I have spent an additional week > rooting out a complex Bulirsch routine from the old Bell lab site and spent > several days relating Jacobi's form to Legendre in order to test against > Abramowitz tables and GSL software for at least real arguments---they now > agree. But Dozier's code does not agree to better than the third or forth > place (I am not saying Dozier's code is wrong yet). > > And I still have no handle on a replacement for Newton-Raphson. What surprises me is that I gave an alternative twice, once in a posting on this list dated 2006-06-13, once in a private email dated 2006-06-23, including a web reference. It's about ACM TOMS algorithm 365 from H. Bach. Fortran code can be found at netlib.org and one can buy a copy of the peer reviewed TOMS article at ACM (acm.org). I am not happy with the method because it is horribly slow, but it is a very robust method, at least if one knows how to tweak its control parameters. As for a replacement of complex Jacobi elliptic functions and a complex elliptic integral: it can be done without too much effort. I provided in the same manner (list posting and private email) references. I do admit that working with complex elliptics is no trivial endeavor without knowledge or experience. _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
On Wednesday 05 July 2006 8:31 am, Oscar van Vlijmen wrote:
> > From: "Gerald I. Evenden" > > > > Stating it another way, my surrender was to the Dozier article that while > > the method sounds interesting the execution falls short. The code had > > several errors and some methodology is in question. I have spent an > > additional week rooting out a complex Bulirsch routine from the old Bell > > lab site and spent several days relating Jacobi's form to Legendre in > > order to test against Abramowitz tables and GSL software for at least > > real arguments---they now agree. But Dozier's code does not agree to > > better than the third or forth place (I am not saying Dozier's code is > > wrong yet). > > > > And I still have no handle on a replacement for Newton-Raphson. > > What surprises me is that I gave an alternative twice, once in a posting on > this list dated 2006-06-13, once in a private email dated 2006-06-23, > including a web reference. > It's about ACM TOMS algorithm 365 from H. Bach. Fortran code can be found > at netlib.org and one can buy a copy of the peer reviewed TOMS article at > ACM (acm.org). > I am not happy with the method because it is horribly slow, but it is a > very robust method, at least if one knows how to tweak its control > parameters. > > As for a replacement of complex Jacobi elliptic functions and a complex > elliptic integral: it can be done without too much effort. I provided in > the same manner (list posting and private email) references. > I do admit that working with complex elliptics is no trivial endeavor > without knowledge or experience. The Bulirsch code seems to work and compares in the real case (x+0i) with the GSL library *and* Abramowich tables--- after a lot of footwork dealing with figuring out how to match up the (&*^(&^ parameters to each. The input parameters were what was really blowing my mind. I ran into another, perhaps the same source, FORTRAN source but found out I did not have a g77 compiler (I had recently changed versions of Linux) and solving that problem introduced more headaches than I wanted at the moment. In my little bench program where I compared GSL, Bulirsch and Dozier I ran into consistant disagreement with Dozier code in the 3rd or 4th decimal place. At the moment, this does not make sense because Dozier's code is at a significant point in Dozier's second phase of the forward solution I do not understand how it could be incorrect and yet yield correct answers in many areas for the overall projection. The following is a short example of a test run: e^2 used: 0.006722670022 x + y*i input: 0.1 0 Dozier e2ci: 0.09999998496533813 0*i Bulirsch elco2 ret: OK elco2: 0.0999988817825081 0*i gsl_sf_ellint_E_e x: 0.09999888178250819 --- err: 4.44104e-17 x + y*i input: 0.2 0 Dozier e2ci: 0.19999988044241 0*i Bulirsch elco2 ret: OK elco2: 0.1999911075211161 0*i gsl_sf_ellint_E_e x: 0.1999911075211167 --- err: 8.88297e-17 x + y*i input: 0.5 0 Dozier e2ci: 0.4999982088519777 0*i Bulirsch elco2 ret: OK elco2: 0.4998667513591754 0*i gsl_sf_ellint_E_e x: 0.4998667513591855 --- err: 2.22222e-16 x + y*i input: 0.3 0.3 Dozier e2ci: 0.3000007837129038 0.2999991577227006*i Bulirsch elco2 ret: OK elco2: 0.3000583008002545 0.2999373655430532*i In the last case, GSL not run because y != 0. That gives you an idea where I now stand with the complex elliptical integral of the second kind (incomplete). The nasty part is that even if this is solved, the problem with the first phase and the root finding looms as an even bigger issue. Sigh. :-( > _______________________________________________ > Proj mailing list > [hidden email] > http://lists.maptools.org/mailman/listinfo/proj -- Jerry and the low-riders: Daisy Mae and Joshua "Cogito cogito ergo cogito sum" Ambrose Bierce, The Devil's Dictionary _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
In reply to this post by OvV_HN
On Wednesday 05 July 2006 8:31 am, Oscar van Vlijmen wrote:
... > What surprises me is that I gave an alternative twice, once in a posting on > this list dated 2006-06-13, once in a private email dated 2006-06-23, > including a web reference. > It's about ACM TOMS algorithm 365 from H. Bach. Fortran code can be found > at netlib.org and one can buy a copy of the peer reviewed TOMS article at > ACM (acm.org). > I am not happy with the method because it is horribly slow, but it is a > very robust method, at least if one knows how to tweak its control > parameters. Not sure what you meant by slow but the example clocked in at 2ms real time and did not show on the millisecond user/system clocks. Sorta looks like a steepest descent routine with no analysis from derivatives so it is likely to be sluggish. So much for Dozier's idea of a speedy method. ;-) I stupidly missed all the control in the file so had to do about three whacks at getting it to compile. Once the junk was removed it went fine. -- Jerry and the low-riders: Daisy Mae and Joshua "Cogito cogito ergo cogito sum" Ambrose Bierce, The Devil's Dictionary _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
Hi, I hope this is not a FAQ or covered in the existing docs; I've
checked and apologize if I missed something obvious. I have a data server which is used to access a number of files of georeferenced data stored in two dimentional arrays. For each array, I know the longitude and latitude values at the corners as well as the projection and datum. Is there a way I can use libproj to map long/lat locations to indices within the arrays? Thanks, James -- James Gallagher jgallagher at opendap.org OPeNDAP, Inc 406.723.8663 _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
James,
What you're trying to do sounds like re-projecting a grid. (A grid just being a two dimensional array of values). I have written code to do this previously. To project a grid you basically have to perform a forward projection of some points around the edge of the input grid and use these to compute the bounding rectangle in the new coordinate system. (Remember some projections will distort the shape of the grid). You then create an output grid in the new coordinate system with the appropriate dimensions. Finally you have to go through each cell in the output grid and use reverse projection to find the corresponding cell in the unprojected grid. Use grid interpolation to pick up the value you need to put in the output grid cell. There's probably already code in the public domain which uses PROJ to do this. Anthony. --- James Gallagher <[hidden email]> wrote: > Hi, I hope this is not a FAQ or covered in the existing docs; I've > checked and apologize if I missed something obvious. > > I have a data server which is used to access a number of files of > georeferenced data stored in two dimentional arrays. For each array, > > I know the longitude and latitude values at the corners as well as > the projection and datum. Is there a way I can use libproj to map > long/lat locations to indices within the arrays? > > Thanks, > James > -- > James Gallagher jgallagher at opendap.org > OPeNDAP, Inc 406.723.8663 > > _______________________________________________ > Proj mailing list > [hidden email] > http://lists.maptools.org/mailman/listinfo/proj > __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
In reply to this post by James Gallagher-4
If your corner coordinates aren't too far apart, then simple linear
interpolation will probably work fine to translate coordinates to array indices. Something like: X_width = Abs( X1 - X0 ) / X_count ; ( X0, Y0 ) minimum x/y Y_width = Abs( Y1 - Y0 ) / Y_count ; ( X1, Y1 ) maximum x/y i = Floor( (Xi - X0) / X_width ) ; ( Xi, Yi ) the coordinate of interest j = Floor( (Yi - Y0) / Y_width ) (Assumes X0 <= Xi < X1 and Y0 <= Yi < Y1) >>> [hidden email] 7/5/2006 1:37:56 PM >>> Hi, I hope this is not a FAQ or covered in the existing docs; I've checked and apologize if I missed something obvious. I have a data server which is used to access a number of files of georeferenced data stored in two dimentional arrays. For each array, I know the longitude and latitude values at the corners as well as the projection and datum. Is there a way I can use libproj to map long/lat locations to indices within the arrays? Thanks, James -- James Gallagher jgallagher at opendap.org OPeNDAP, Inc 406.723.8663 _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
In reply to this post by Anthony Dunk
On Jul 5, 2006, at 9:24 PM, Anthony Dunk wrote: > James, > > What you're trying to do sounds like re-projecting a grid. (A grid > just > being a two dimensional array of values). I have written code to do > this previously. I think what I need is simpler, but maybe not. What I'm doing is writing a function that can be run by a data server that will return a region of interest. The ROI, suplied by a client, will specified by long/lat upper-left and lower-right corner points. What I need to get are the actual array indexes that correspond to those long/lat locations. For some of the data, the coverage is the 'global' (which means 80N to 80S and 180W to 180E) while for the rest it a smaller region. Metadata associated with these data files are spotty and often somewhat odd. For example, one data set lists longitude as running from 20 to 380 degrees. Before I look more at libproj, I should ask if it will handle such oddities? I assume so, but ... Thanks, James > > To project a grid you basically have to perform a forward > projection of > some points around the edge of the input grid and use these to compute > the bounding rectangle in the new coordinate system. (Remember some > projections will distort the shape of the grid). > > You then create an output grid in the new coordinate system with the > appropriate dimensions. > > Finally you have to go through each cell in the output grid and use > reverse projection to find the corresponding cell in the unprojected > grid. Use grid interpolation to pick up the value you need to put in > the output grid cell. > > There's probably already code in the public domain which uses PROJ to > do this. > > Anthony. > > --- James Gallagher <[hidden email]> wrote: > >> Hi, I hope this is not a FAQ or covered in the existing docs; I've >> checked and apologize if I missed something obvious. >> >> I have a data server which is used to access a number of files of >> georeferenced data stored in two dimentional arrays. For each array, >> >> I know the longitude and latitude values at the corners as well as >> the projection and datum. Is there a way I can use libproj to map >> long/lat locations to indices within the arrays? >> >> Thanks, >> James >> -- >> James Gallagher jgallagher at opendap.org >> OPeNDAP, Inc 406.723.8663 >> >> _______________________________________________ >> Proj mailing list >> [hidden email] >> http://lists.maptools.org/mailman/listinfo/proj >> > > > > > __________________________________________________ > Do You Yahoo!? > Tired of spam? Yahoo! Mail has the best spam protection around > http://mail.yahoo.com > _______________________________________________ > Proj mailing list > [hidden email] > http://lists.maptools.org/mailman/listinfo/proj -- James Gallagher jgallagher at opendap.org OPeNDAP, Inc 406.723.8663 _______________________________________________ Proj mailing list [hidden email] http://lists.maptools.org/mailman/listinfo/proj |
Free forum by Nabble | Edit this page |