Bug #831

-D in grdproject

Added by Mark over 1 year ago. Updated 9 months ago.

Status:FeedbackStart date:2015-12-05
Priority:NormalDue date:
Assignee:Paul% Done:


Target version:Candidate for next bugfix release
Affected version:5.2-svn Platform:Mac OS X


I am using grdproject for a forward projection. The following command works fine, setting the output grid to the same dimensions as the input grid.

gmt grdproject in.grd -Gout.grd -Jt-127.5/44.25/1:1 -Fk -C -V

However, the command fails when I set -D to get a denser grid, in this case a 1 km grid spacing:
gmt grdproject in.grd -Gout.grd -Jt-127.5/44.25/1:1 -Fk -C -D1k -V

The explanation of the -D option indicates that it is available to set the grid spacing of the output grid, but the examples suggest that it only works for inverse projections. If this is the case, then it should be stated in the -D option.

In any case, it my sense that there is no way to specify the grid spacing for a forward projection.


#1 Updated by Paul over 1 year ago

  • Status changed from New to Feedback
  • Assignee set to Paul
  • Target version set to Candidate for next bugfix release

Think there are two problems here that need to be addressed:

  1. The -D option does not take units since you are projecting from degrees to km. Of course, the man page says otherwise since it is using a shared documentation template. This confusion we can fix. So in your case you only want to say -D1 (since that means 1 km as your chosen using via -Fk).
  2. However, grdproject still honors the projected min/xmax/ymin/ymax and then adjusts dx,dy to fit, given nx,ny and that region. I just made a dummy grid since yours were not attached but my -D1 ended up becoming dx = 1.00019738 and dy = 1.00003761 to keep the region the same. Presumably, when you use -D you wish to round the projected region to the nearest multiple of that increment in the outward sense, so that the rounded region is >= the original region. I cannot see a situation when we want to keep the original region and readjust the stated increment (as it does now), but maybe the best option is to provide a way to do either. Perhaps adding new option -A to mean: In conjunction with -D, round projected region to a whole multiple of the chosen increments?

Just so on same page: I made a test grid like this: gmt grdmath -R-140/-120/40/50 -I2m X = in.grd and projected using your arguments to get

out.grd: x_min: -1066.17081203 x_max: 639.16571785 x_inc: 2.84222755 name: km nx: 601
out.grd: y_min: -472.57889775 y_max: 714.46574895 y_inc: 3.95681549 name: km ny: 301

and using -D1 only changes nx,ny and adjusts dx,dy, but keeps the same domain:

out.grd: x_min: -1066.17081203 x_max: 639.16571785 x_inc: 1.00019738 name: km nx: 1706
out.grd: y_min: -472.57889775 y_max: 714.46574895 y_inc: 1.00003761 name: km ny: 1188

With the proposed -A -D1 I would expect

out.grd: x_min: -1067 x_max: 640 x_inc: 1 name: km nx: 1708
out.grd: y_min: -473 y_max: 715 y_inc: 1 name: km ny: 1189

Would this be what you would prefer?

#2 Updated by Mark over 1 year ago

Got it... -D picks up its units from -Fk. Some details in the documentation would solve all my interests.
Thanks a bunch...

#3 Updated by Paul over 1 year ago

OK, but would you not also want to adjust -R to a multiple of -D units?

#4 Updated by Mark over 1 year ago

Yes, that would be good, and consistent with other gmt commands...

#5 Updated by Joaquim 9 months ago

But because -D, like -I, calls gmt_getinc() one can do -D1= and get

out.grd: x_min: -1068.57902517 x_max: 640.420974829 x_inc: 1 name: km n_columns: 1710
out.grd: y_min: -493.352580242 y_max: 692.647419758 y_inc: 1 name: km n_rows: 1187

Also available in: Atom PDF