Mapproject using forward projection with oblique Mercator
|Target version:||Candidate for next bugfix release|
|Affected version:||5.1.x||Platform:||Mac OS X|
I have been using mapproject with the oblique Mercator projection for some plate motion calculations. I ran into two problems that are either bugs or unexplained features.
The problems are demonstrated with some simple examples in the attached shell script. Of course, these issues may be due to an incorrect vision of how mapproject should work with an oblique Mercator projection. In particular, I am using -Joc followed by a lon,lat defining the projection center, and lon,lat defining the oblique projection pole. My expectation is that p,q would be defined in a right-handed sense relative to the projection center and projection pole, and p,q would be equal to 0,0 at the projection center. What I find in the more complicated examples I am working, and demonstrated in the simple examples attached here, is that the p,q coordinates are the same for an oblique projection pole and its antipodal pair, and that the origin for the p,q coordinate frame is always fixed at the intersection between oblique equator and the oblique central meridian, regardless of the oblique latitude of the specified projection center.
I can work my way around these issue, but I suspect that they are due to oversights in the calculations. It took me a long time to sort all this out, so it might help to either change the code or to provide a bit more explanation in the documentation.
Make sure mapproject -C gives (0,0) when origin is given for oblique Mercator, issue #829
#1 Updated by Paul almost 2 years ago
- Status changed from New to Feedback
- Assignee set to Paul
- Target version set to Candidate for next bugfix release
Hi Mark. So mapproject was never written with the intention to serve plate tectonics (but of course aspects of it can). Perhaps you could state your actual need and maybe there are other tools that can better serve those needs, e.g., gmtvector or the spotter supplement.
#2 Updated by Mark almost 2 years ago
The plate tectonics application is not the point. My point is that given the specification of the projection center and the projection pole, one should be able to anticipate the outcome for p and q. Right now, the documentation does not provide this information. I would expect that the approach would follow some common reference, such the USGS Synder approach.
#3 Updated by Paul almost 2 years ago
Will see what can be done, but be aware AGU is looming... In the last release we did add uppercase A,B,C, which means Upper-case A|B|C removes enforcement of a northern hemisphere pole. So I tried your cases with -JoC but it did not help for these cases - it may only be effective for plots.
#5 Updated by Paul about 1 year ago
- Status changed from Feedback to In Progress
- % Done changed from 0 to 50
Apparently AGU is now over and time to look at this again...
So baby-steps: I agree that mapproject -C should yield (0,0) when the projection origin is given. I consider that a bug and it has now been fixed in r16819. This only affects mapproject -C output since plots don't care.
As for having p go negative with increasing longitude: This is not so clear. It is certainly true and the case for project but mapproject is a bit different. I will do some experiments to see if this could behave more sensibly or at least update the documentation on that part.
#6 Updated by Paul about 1 year ago
- Status changed from In Progress to Resolved
- % Done changed from 50 to 100
More baby steps: r16822. Now p,q behaves the way you expect it: (p,q) = (0,0) at the origin, q increases toward oblique pole, and p increases towards the right when looking toward the pole. Let us know fi there are any issues.