How to use grdproject with oblique mercator?

Hi all,

I’m trying to use the oblique mercator projection to get an equal x,y grid in meters over the Gakkel Ridge using the IBCAO arctic bathy grid. Trying to center the projection at 8,85 with the equator on an azimuth of 46. grdproject runs, but gives me this error:
grdproject [WARNING]: x_inc does not divide 180; geographic boundary condition changed to natural.
which I thought the +e on the -D option would take care of.

In any case, it creates a grid, but doesn’t add the Z values. Also, the grid dimensions are far too large for the area I selected.

test1.grd: Title: Produced by grdproject
test1.grd: Command: grdproject IBCAO_V3_30arcsec_RR.grd -R-12/30/83/86.5 -Joa8/85/46/1:10000 -Fe -D3000+e -Gtest1.grd
test1.grd: Remark:
test1.grd: Gridline node registration used [Cartesian grid]
test1.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
test1.grd: x_min: 0 x_max: 4671000 x_inc: 3000 name: m n_columns: 1558
test1.grd: y_min: 0 y_max: 4422000 y_inc: 3000 name: m n_rows: 1475
test1.grd: z_min: 0 z_max: 0 name: z
test1.grd: scale_factor: 1 add_offset: 0
test1.grd: format: netCDF-4 chunk_size: 130,135 shuffle: on deflation_level: 3

Am I trying to use this wrong, or misunderstanding the man page?

Rick Hagen
Naval Research Lab
Washington, DC

OK, well to reply to myself, part of my problem is that I mis-specified -R, which should be -R-12/83/30/86.5+r

My entire command is now: grdproject IBCAO_V3_30arcsec_RR.grd -R-12/83/30/86.5+r -Joa8/85/46/1:10000 -Fe -D3000+e -Gtest1.grd

I now get z-values in my grid, and it is in the right area. However, I expected that by specifying the oblique equator azimuth along the ridge the grid would be rotated so that the ridge would be horizontal when plotted. That’s not happening, so I guess I mis-understood.

Well, to reply to myself again: I am baffled. I am able to reproduce the oblique mercator plot of New Zealand from the GMT cookbook, but can’t get my area to plot correctly using essentially the same projection parameters.

The cookbook script:
pscoast -R-1000/1000/-500/500+uk -JOa173:17:02E/41:16:15S/35/3i -G111/111/111 -B5g5 -P >

My script:
grdcontour IBCAO_V3_30arcsec_RR.grd -R-200/200/-200/200+uk -JOa8/85/46/4i -C100 -B5g5 -P >

The plot of my area is centered correctly, but the 46 degree azimuth of the equator is not rotated to the horizontal.

Help Please!

I don’t know how to fix this. I just notice that the problem seems to be related to the proximity of the pole. If you plot a little to the south it looks ok:

gmt basemap -R-200/200/-200/200+uk -JOa8/80/46/4i -B5g5 -T8/80/3c -png test

Yes Marcelo, just before you posted I ran a test at a random spot near 70N and it worked. Glad to have independent confirmation. Maybe a limitation of the projection? or a bug?

Just to update this. Turns out there is a bug with -Joa, it doesn’t work north of 83N. I was able to get this to work using -Job:
grdproject IBCAO_V3_30arcsec_RR.grd -R-12/83/30/86.5+r -Job8/85/15.1/85.5/1:1 -Fe -D3000+e -Gtest1.grd