I use mapproject to calculate coordinates for lower left and upper right corner for plotting a rectangular map similar to gmt Example 28 (some of the plotted data are projected according to EPSG:3035, other are geographical lon/lat degrees). When I specify coordinate system using EPSG code, mapproject produces output as expected, but also displays an error message 3 times:
echo 4200000 3360000 | gmt mapproject -JEPSG:3035/1:1 -I -F
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
8.18360619615 53.3343880853
If I specify EPSG:3035 coordinate system manually, there is no error message:
echo 4200000 3360000 | gmt mapproject -Ja10/52/1:1 -F -C4321000/3210000 -I -R-16.1/32.88/40.18/84.17
8.1836061961 53.3343880946
There is also a very small difference in the latitude, not critical to my application.
Is this anything of importance? The repeated error message is quite annoying.
gmt 6.4.0, miniconda environment on Linux.
The difference is that when you do the proj specification it runs via GDAL and those ERROR 1 messages are from inside GDAL, not GMT, and outside our control. Maybe @Joaquim has some more info on this.
Thanks @pwessel for clarifications.
Below is an extract from debug messages. There is a strange line mapproject [DEBUG]: Projected values in meters: -9.00996e+06 9.00996e+06 -9.00996e+06 9.00996e+06 right before the triplicate error message. I have no idea where those values come from. As well as mapproject [INFORMATION]: Transform 0/0/-90/90 <- -9009964.7612/9009964.7612/-9009964.7612/9009964.7612 [m] that is right after the errors. I haven’t specified these particulate values on my input.
gmt [DEBUG]: Revised options: -J+EPSG:3035/1:1 -I -F -Vd
mapproject [DEBUG]: History: Process -J+EPSG:3035/1:1
mapproject [INFORMATION]: Converted to -J syntax = -Ja10/52/1:1
mapproject [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
mapproject [INFORMATION]: Processing input table data
mapproject [DEBUG]: Reset MAP_ANNOT_OBLIQUE to anywhere
mapproject [DEBUG]: Projected values in meters: -9.00996e+06 9.00996e+06 -9.00996e+06 9.00996e+06
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
mapproject [INFORMATION]: Transform 0/0/-90/90 <- -9009964.7612/9009964.7612/-9009964.7612/9009964.7612 [m]
mapproject [DEBUG]: gmtapi_init_import: Passed family = Data Table and geometry = Point
mapproject [DEBUG]: gmtapi_init_import: Added 0 new sources
mapproject [DEBUG]: Object ID 0 : Registered Data Table Stream 7f616ce7d980 as an Input resource with geometry Point [n_objects = 1]
mapproject [DEBUG]: gmtapi_init_import: Added stdin to registered sources
No idea why. Same numbers when doing the direct projection (note that you don’t need the 1:1 -F)
echo 8.18360619615 53.3343880853 | gmt mapproject -JEPSG:3035 -Vd
...
mapproject [DEBUG]: Reset MAP_ANNOT_OBLIQUE to anywhere
mapproject [DEBUG]: Projected values in meters: -9.00996e+06 9.00996e+06 -9.00996e+06 9.00996e+06
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
mapproject [INFORMATION]: Transform 0/0/-90/90 -> -9009964.7612/9009964.7612/-9009964.7612/9009964.7612 [m]
I guess it’s doing some other operation inherited from the GMT proj machinery and in this case it triggers those GDAL errors. For example this case that converts to UTM29 doesn’t print the errors but also display that type of info
echo -8 40 | mapproject -J32629 -Vd
mapproject [INFORMATION]: Converted to -J syntax = -Ju29/1:1
mapproject [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
mapproject [INFORMATION]: Processing input table data
mapproject [DEBUG]: Reset MAP_ANNOT_OBLIQUE to anywhere
mapproject [DEBUG]: Projected values in meters: 1.50549e+06 1.61848e+06 0 112247
mapproject [INFORMATION]: Transform 0/1/0/1 -> 1505491.99671/1618480.92034/0/112246.665666 [m]
The Converted to -J syntax step should not be needed in mapproject but it is in grdimage for example because we need that info to create the map frames for the different projections. This is the reason why we cannot make figures in all the PROJ supported projections. Well, in fact we can, but we can’t create the frames for the GMT unknown ones.
Thanks @Joaquim for the explanations.