@Joaquim and @pwessel, thanks for your replies. Your discussion really helped clarify how GMT deals with coordinate systems. It is counter-intuitive to me as it seems to use a different concept from others QGIS or the recent cartopy…
In the ex28, it says
“The mistake many GMT rookies make is to specify the UTM projection with their UTM data. However, that data have already been projected and is now in linear meters. The only sensible way to plot such data is with a linear projection, yielding a UTM map.”
I think I’m that rookie
To double check my understanding, here is my summary. The plotting functions in GMT needs geographic coordinates (lon/lat) by default all the time (either table data or grid data). The -J in GMT functions is to tell these functions how to put given lon/lat onto a projected planar surface, but NOT to tell these functions which projection system your input coordinates are based on. When your input coordinates are already in projected system, you simply use a planar coordinate system (linear projection) to plot them like you would plot any two-dimensional data points.
Since geographic coordinates (lon/lat) also depends on which ellipsoids they are based upon, it is the users’ responsibility to supply lon/lat in the same ellipsoid as the projection system specified by -J.
The GMT does not check the coordinate system of the input table data file, it will simply assume x/y in a table data file are lon/lat and try to plot them in the specified projection by -J. I confirmed this by trying to plot data from a OGR/GMT ASCII file in the UTM system onto the projection of its UTM zone (I convert a GPKG file in UTM to GMT file using
gmt plot fails. But the GMT does recognize coordinate systems in the grid data files. I confirmed this by trying to plot a Geotiff image in the UTM system onto the projection of its UTM zone.
gmt grdimage does check the coordinate system of an input grid data file and do not simply assume lon/lat coordinates of input grids. However, it won’t plot the image boundary correctly. However, when I use
gdalwarp to reproject the image into a geographic system (EPSG:4326, WGS84) and then
gmt grdimage will plot the image boundary correctly. I think the projection/reprojection procedure has some precision issues? I’ll follow up with the figure plots I generated from the above two cases. Gotta go now.