Grdimage of projected grid fails

Q&A
1

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 tmp.nc`

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

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

2

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.

3

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.

4

Yes, please!

5

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

2 Likes
6

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!

7

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.