My result is: UTM data works, Web Mercator produces two separate maps with an offset, EPSG:3057 (ISN93 / Lambert 1993 - EPSG:3057) produces map but geographical coordinates are wrong, like the projections fast easting and northing is added to the geographical coordinates.
is this expected behaviour, that is, only UTM data can be overlaid by geographical data like in ex28?
my minimal test script and plots below. Would be grateful for explanations - whether this should work or not or I am doing something radically wrong.
GMT version:gmt-6.5.0.dev8+5496e05-h02197e7_0.conda in a mamba environment, Linux.
More like unimplemented features. We never used -JEPSG etc until @joaquim started to play with it so you might have more success with using the equivalent GMT projection syntax. At any rate, need @Joaquim’s input on this.
Unfortunately, for Web Mercator projection there is no equivalent GMT projection syntax.
For EPSG:3057 there is one, -Jl... for Lambert conic conformal. It also works if I define -Jl-19/65/64.25/65.75 using a proj4 string '+proj=lcc +lat_0=65 +lon_0=-19 +lat_1=64.25 +lat_2=65.75'. proj4 equivalent of the EPSG:3057 did not work, like x and y offsets somehow not taken into account.
Web Mercator is just a Mercator but one that uses a sphere instead of ellipsoid. It can be reproduced in GMT by setting the ellipsoid to sphere but need ofc to use the same radius as used in the corresponding EPSG code.
This whole think look more a modern mode issue. Using either EPSG’s or GMT -J syntax gives the close thing.
(Edited the code below, used wrong coordinates in the first version of this comment)
Great, thanks for showing the corresponding GMT syntax for EPSG:3857!
How about EPSG:3057? Please see below.
$ echo 239980 309980 | gmt mapproject -JEPSG:3057/1:1 -C -F -I -V
mapproject [INFORMATION]: Converted to -J syntax = -Jl-19/65/64.25/65.75/1:1
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Central meridian outside region
mapproject [INFORMATION]: Transform 0/1/0/1 <- 3385496.72497/3614199.79014/-8113992.43392/-7886080.5878 [m]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
-24.1704616471 63.2001617189
...
$ echo 239980 309980 | gmt mapproject -Jl-19/65/64.25/65.75/1:1 -C -F -I
-13.3270679583 67.6868858416
$ echo 239980 309980 | awk '{print $1-500000, $2-500000}' | gmt mapproject -Jl-19/65/64.25/65.75/1:1 -C -F -I
-24.1785832117 63.1969909697
Thanks for the explanations re. mapproject. It seems like my coordinate transformations using mapproject were correct after all.
I think my problem is with plotting. I remember you explained me before that using EPSG definitions or PROJ strings are only partly supported for the map frame drawing, so this is probably it.
I use an -R region, but another region appears on the map if I specify -J as EPSG. I use -JEPSG:3057/1:1280000 -R-24.1704616471/63.2001617189/-13.1504259197/66.5062833504+r, but the map is plotted like another region has been specified, something like -R-31.../-58.../-23.../62...+r
The thing which drove me crazy is that I have to use -JEPSG:3057 for the coord transformations but -Jl-19/65/64.25/65.75 for plotting. I don’t understand how to do the trick with EPSG:3857 (web mercator)
By “the trick” I probably mean is this if I understand it correctly:
Is there a way to force ellipsoidal (or other) equations when doing actual plotting with some actual region of interest?
I think the longitude range in -R must be < 10 degrees to force ellipsoidal equations.
For this type of maps it is very confusing (as you found) to mix GMT and PROJ syntaxes. The problem with the frames is that for one side, it is necessary to exactly reproduce the GMT -J syntax from PROJ or EPSG and next that projection must exist in GMT. The only sane way is to find out manually the EPSG/PROJ → -J correspondence and use only it in the plotting commands.
Note that for mapproject and EPSG nothing of the above bothers because projection calculations are done in PROJ without any need for a GMT conversion.
