Error when trying to use projected units in -R by +ue along with EPSG code for -J

Hi, I’m a new to GMT and feeling like on a steep learning with all the syntax of GMT 6.1.0. Here I’m trying to use projected units (meter) in -R by +ue along with EPSG code for -J (UTM zone 33). But I can’t make it work.

Simply using lat/lon in degree works,
gmt basemap -JEPSG:32633+width=16c -R12:48:41.52/53:41:33.11/13:32:14.72/54:8:37.29+r -B -B+glightred+t"My first plot"

But the following won’t work …
gmt basemap -JEPSG:32633+width=16c -R355502.500/5951537.500/404462.500/6000497.500+r+ue -B -B+glightred+t"My first plot"

I’ve looked through the documentation on -R but I can’t see what went wrong.

Thanks,

Zhan.

I do not think the proj6 way of specifying -J works for plotting. I think you need to use -JU with the proper UTM zone. Also you definitively need to remote the +r since you are giving min/max here, not corners.

Hi @pwessel, thanks for your help here. What if I wanna use another projection, EPSG:3035, LAEA Europe. It’s a Lambert_Azimuthal_Equal_Area projection. I know I can potentially configure it with -JA but it’d be much more convenient and less prone to typos to use EPSG code … The current GMT doc page says it supports EPSG code. So where can I use the EPSG code then?

Also, I double checked, the coordinates are in pairs of LL and UR corners :slight_smile:

Zhan.

OK, I was wrong about your +r. We need @Joaquim to comment on this since he did the EPSG stuff.

Hmm, this not about EPSGs. In GMT when one want to do a map projection we give coordinates in geographics and specify the projection (either with classic GMT syntax, PROJ4 syntax or EPSG).
When data is already projected, like in your second example, than we simple use a linear projection, normally -JX=16c

That’s when I think the GMT in Julia makes life a lot easier.

True if the data is in utm meters, but you can still specify the region with utm coordinates provided you append +ue (we compute the corresponding lon,lat under the hood.

Ah, yes ex28. Then I guess it doesn’t know about the EPSG setting at the time of that under the hood operation.

@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 :sweat_smile:

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 gdal_translate). 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.

Thanks again,

Zhan.

Good questions and we should have a kind of FAQ to address them.

(first paragraph) Basically yes.

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.

Yes, but … you can also tell the plotting functions what ellipsoid to use. See PROJ_ELLIPSOID

Third paragraph touches several points. GMT does not generally knows the coordinate system of tabular data. The particular .gmt format may have that info but I don’t think we use it anywhere. So basically, no, you cannot do re-projections on the fly, must use mapproject.
Grids may carry projection info but again this info is not used in a generalized way. This is a part where we clearly need to improve and the basic machinery is already there. For example the Julia version I pointed you is able to read OGR formats directly but not the base GMT.

Another thing that you must be warned is that not all mapping projections (for map plotting) work well with EPSG or the PROJ4 syntax. The arrays projection and image creation do but the image frames may be fail to plot correctly. For example you can make plots with unsupported projections in GMT (but supported in PROJ) but frames will show up funny.

Aha, PROJ_ELLIPSOID! There are many more to learn about GMT. But I haven’t found ways to designate the datum upon which lon/lat are defined in the plotting function. Only mapproject has an option -E to designate datum.

Another thing that you must be warned is that not all mapping projections (for map plotting) work well with EPSG or the PROJ4 syntax. The arrays projection and image creation do but the image frames may be fail to plot correctly. For example you can make plots with unsupported projections in GMT (but supported in PROJ) but frames will show up funny .

Ah now I see why the image frames I plotted look funny. It is the EPSG or the PROJ4 syntax. If I want some “perfection” for image frames, I need to use -J with detailed projection parameters according to the projection type rather than a simple EPSG code or PROJ4 string.

Thanks again for your explanation! Indeed more improvements are needed. But even given current functionality, I’ve found GMT quite helpful being a pure command-line mapping tool to batch generate maps for productivity and reproducibility!

Zhan.

You can use gmt set PROJ_ELLIPSOID name or --PROJ_ELLIPSOID=name in any of the map plotting programs

Hmmm,but PROJ_ELLIPSOID only sets ellipsoid, doesn’t ? It does not seem to set up datum though to me. Sometimes, I work with data in either ETRS89 or NAD83 datum, which both base upon the ellipsoid GRS80. PROJ_ELLIPSOID has the option GRS80 but not any specific datum that is based on this ellipsoid. Maybe I’m missing some syntax to designate a datum?

Yes, that is true and then the user must take of the datum transforms by himself. But if you are interested in mapping at such large scales where that matters then mind you that even WGS84 is not one but 6 different realizations. If you are interested, proj added time and deformation models in its transforms. See https://github.com/OSGeo/PROJ/pull/2206

Thanks for pointing out this new function in PROJ!