gmt icon indicating copy to clipboard operation
gmt copied to clipboard

-Jepsg where epsg is an UTM with negative longitudes fails in classic but only for .ps in modern.

Open joa-quim opened this issue 2 years ago • 3 comments

The modern case first. This fails. No PS error but figure is blank

gmt coast  -R-81.8/-66.9/-1.0/15.0 -J32618 -Baf -BWSen -W0.5 -Da -ps lixo

but this, or any other format, works.

gmt coast  -R-81.8/-66.9/-1.0/15.0 -J32618 -Baf -BWSen -W0.5 -Da -eps lixo

Now the classic mode case. Since here we only have output in .ps it always fails. But using the -J GMT syntax, works. That is, this works. (epsg:32618 <==> JU18)

gmt pscoast  -R-81.8/-66.9/-1.0/15.0 -JU18/14c -Baf -BWSen -W0.5 -Da -P > lixo.ps

Several debugging hours later I found why. The gmtproj_utm() function (called by gmt_geo_to_xy) has

if (!GMT->current.proj.north_pole) (*y) += GMT_FALSE_NORTHING;	/* For S hemisphere, add 10^7 m */

but the corresponding GDAL function that is used for the UTM transformation doesn't have that offset and thus the computed ps units are negative and points fall off cliff.

So, one of the problems is identified though I don't know a cure since the transformation for the epsg case takes place inside the GDAL lib.

The second problem is: HTF the modern mode is insensitive to problem one when not using the .ps format????

joa-quim avatar Aug 31 '23 22:08 joa-quim

I could not easily find where the epsg:32618 -> JU18 translation happens. But I would think a solution like this should work

  1. Add a new parameter bool GMT->current.proj.north_offset.
  2. When UTM, and north pole, set it to true unless it was the proj translation
    1. Replace the 4 places in gmt_proj.c where GMT_FALSE_NORTHING appears.

PaulWessel avatar Sep 01 '23 07:09 PaulWessel

I could not easily find where the epsg:32618 -> JU18 translation happens.

gmtinit_parse_proj4 calls gmt_importproj4

And how would your recipe be applied? Remember that for mapping we need the equivalence between proj4 <-> GMT because of the map frames. These can only be plotted from the GMT projection machinery but the inners of the maps are computed from proj. This sad fact is what prevent us to make maps in any of proj supported projections.

joa-quim avatar Sep 01 '23 11:09 joa-quim

Since modern mode has a 5 inch buffer in that 10x10 m paper it adjusts up to 5" based on negative coordinates. Try adding a larger size (e.g/18c) and it gets clipped at the base.

Seems it ought to be possible to set a flag that means GDAL UTM is used and somewhere in geo_to_xy we insert the if-test to add to y?

PaulWessel avatar Nov 24 '23 15:11 PaulWessel