Grdimage of projected grid fails

I have a raster that is already in projected units (meters) that I would like to plot using grdimage by using the +ue suffix to the -B string.The following code fragment:

bounds=`gmt grdinfo -Ir`

gmt psbasemap $bounds+ue -JA-100/50/10 -B0 -P -K > $output

gmt grdimage -R -J -C$cpt -O -P -K >> $output

Fails with the error

grdimage (gmtapi_import_grid): Subset x range <= 0.0

The header of the grd file looks ok:

gmt grdinfo Title: Produced by grdmath Command: gmt grdmath ../ 100 MUL = Remark: Pixel node registration used [Cartesian grid] Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7 x_min: -2050000 x_max: -155000 x_inc: 5000 name: longitude [meters] n_columns: 379 y_min: -2035000 y_max: 550000 y_inc: 5000 name: latitude [meters] n_rows: 517 v_min: -1.19209289551e-05 v_max: 99.9270706177 name: Barren land fraction [fraction] scale_factor: 1 add_offset: 0 format: netCDF-4 chunk_size: 190,130 shuffle: on deflation_level: 3 Default CPT: 

Any idea if I am misunderstanding the use of the projected grid coordinates here?

The workaround I have at the moment is to plot the raster with grdimage using a cartesian grid (-JX) and then calling psbasemap again with -R+ue so I can plot vectors that are in lon-lat. This works, but requires calculating the aspect ratio of the plot in advance in order to pass the correct plot size to -JX. Is there any other way to tell grdimage that a grd file is already in the projected coordinate system and so not to reproject?

Thanks, Jed

Hi, can you try to « reproject » it anyway using grdproject.
That would convert the meters to degrees from a reference point. It’s one line and no calculations are involved on your end.

Normally, you would plot your projected grid with -Jx1:xxxxxxxx and then overlay a basemap using -Ja-100/50/1:xxxxxx using the -Rxmin/ymin/xmax/ymax+r region and whatever -B you want.

I notice that it would be nice if grdinfo could return the xmin/ymin/xmax/ymax+r string when requested so I think I will add -Io to return the lower left/upper right +r values.

Yes, please!

The -Io has now been implemented, tested, and merged into master and can thus be used if you build from git source.


Ok I see, yes using the 1:xxxx instead of absolute plot size avoids the need to calculate the plot aspect ratio. Using projected coordinates as -Rxmin/ymin/xmax/ymax+ue is in my mind preferable to specifying the lower left/upper right plot +r coordinates in decimal degrees but either way it works.

Thanks for the replies!

You still use projected coordinates but the -Io will format that info the way it needs to be passed (with +r). No decimal degrees involved.