Overlaying geographical data over projected cartesian data: UTM works, Web mercator doesn't, EPSG:3057 doesn't

Hi, I am trying to use overlaying of geographical data over projected Cartesian data for some of my map plotting, the same way as in gmt example 28: (28) Mixing UTM and geographic data sets — GMT 6.4.0 documentation http://docs.generic-mapping-tools.org/4.5/gmt/html/GMT_Docs.html#x1-1630007.28

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.

# UTM zone 5

## EPSG:3057

## Web Mercator

X_MIN=$(echo $ROI | awk 'BEGIN {FS="/"} {print $1}')
Y_MIN=$(echo $ROI | awk 'BEGIN {FS="/"} {print $3}')
X_MAX=$(echo $ROI | awk 'BEGIN {FS="/"} {print $2}')
Y_MAX=$(echo $ROI | awk 'BEGIN {FS="/"} {print $4}')
echo "ROI is: $ROI"

LL=$(echo "$X_MIN $Y_MIN" | gmt mapproject -J"$PROJ"/1:1 -F -I --OUTPUT_DEGREE_FORMAT=D | awk '{printf "%s/%s\n", $1, $2}')
UR=$(echo "$X_MAX $Y_MAX" | gmt mapproject -J"$PROJ"/1:1 -F -I --OUTPUT_DEGREE_FORMAT=D | awk '{printf "%s/%s\n", $1, $2}')
printf "Geographical ROI is: -R${LL}/${UR}+r\n"

gmt begin test_frame_PROJ_AFTER_CARTESIAN png -C
  gmt basemap -Bxag -Byag -BSW -Jx$SCALE -R$ROI \
                  --MAP_FRAME_TYPE=plain --MAP_GRID_CROSS_SIZE_PRIMARY=0.25
  gmt basemap -BNE -Bxag -Byag -J"$PROJ"/$SCALE -R${LL}/${UR}+r \
gmt end

UTM grids get overlaid as expected:

EPSG:3057 with “wrong” geographical region, geographical coordinates off by exactly the false northing and easting of EPGS:3057:

Web Mercator with big vertical offset between two maps:

This looks like a bug to me. Do you mind opening a bug report on GitHub?

Hi and thanks for your response.
I don’t mind but unfortunately I don’t know how to do things on GitHub.

OK. Let me ping @pwessel first to confirm if it’s a bug.

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.

And big thanks to @joaquim for enabling this at least partially.

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.

echo 239980 309980 | gmt mapproject -J3054 -I -V
mapproject [INFORMATION]: Converted to -J syntax = -Ju26/1:1
-29.3387226553  2.80129334688
echo 239980 309980 | gmt mapproject -Ju26/1:1 -C -F -I
-29.3388616002  2.80211558448

The difference is for sure due to the fact that we are not doing exactly the same thing, as EPSG 3054 is in fact

julia> epsg2proj(3054)
"+proj=utm +zone=26 +ellps=intl +towgs84=-73,47,-83,0,0,0,0 +units=m +no_defs"

But all this to say that the large differences cannot be coming from the mapproject outputs.

Hmmm, I made a EPSG number confusion. Here is the comparison for the Web Mercator

C:\v>echo 239980 309980 | gmt mapproject -J3857 -I
2.15577701883   2.78350216186

C:\v>echo 239980 309980 | gmt mapproject -Jm0/0/1:1 -C -F -I --PROJ_ELLIPSOID=6378137 -R0/3/0/3
2.15577701883   2.78350216186

(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

PROJ4 definition of EPSG:3057:

+proj=lcc +lat_0=65 +lon_0=-19 +lat_1=64.25 +lat_2=65.75 +x_0=500000 +y_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs
julia> epsg2proj(3057)
"+proj=lcc +lat_0=65 +lon_0=-19 +lat_1=64.25 +lat_2=65.75 +x_0=500000 +y_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
echo -7088277.75767 5567888.65785 | gmt mapproject -I -J3057
-135.259272902  24.0368999372

Now for the -J we need to specify a narrow -R to force the ellipsoidal equations. I hate this and thought it was relaxed in GMT6.5 (@pwessel ?)

echo -7088277.75767 5567888.65785 | gmt mapproject -C500000/500000 -F -I -Jl-19/65/64.25/65.75/1:1 -R-136/-134/0/30
-135.259272902  24.036899937

No need for the /1:1 -C -F when using EPSG or +proj=...

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

gmt basemap -JEPSG:3057/1:1280000 -R-24.1704616471/63.2001617189/-13.1504259197/66.5062833504+r -BWSne -Ba1g1 -png quick-3057

If I specify -Jl... manually, the correct region appears on the map:

gmt basemap -Jl-19/65/64.25/65.75/1:1280000 -R-24.1704616471/63.2001617189/-13.1504259197/66.5062833504+r -BWSne -Ba1g1 -png quick-jl

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 showed above how to do for Web Mercator

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.

Thanks again! Explicitly configuring the projection ellipsoid using --PROJ_ELLIPSOID=6378137 seems fixed basemap frame creation for EPSG:3857 (web mercator)!