How to generate a lat/long geotiff with grdimage

No need for the gdal_edit step with the psconvert way.

Good that it works now. Indeed it’s strange why it failed before. Worked for me on WSL too.

Glad that it works now

If I don’t apply the gdal_edit step, the geotiff crs is unknown (from gdalinfo), but the coords are correct. After that step, gdalinfo gives the correct crs data, so I prefer to apply that step for a complete solution.

Cheers

And not as rosy as I thought… and I think this identifies the problem (or helps at least)…

I’m using on Linux (as suggested by Joaquim):

gmt grdimage @earth_gebco_05m_g -P -J4326 > world.ps
gmt psconvert world.ps 1W+g
# add the edit to set CRS
gdal_edit.py -a_crs EPSG:4326 world.tiff

It works with no errors for resolutions lower than or equal to _06m_

This is a single grdfile download. I start getting my errors at higher resolutions than this.

At ~_05m_g it downloads multiple grids to achieve the global coverage & I get the errors described above. The problem seems to be in the merging of grids, wherever that occurs.

I get this same error under Linux & Windows.

Could you please try this and see if it still works for you…

Brent

I did exactly like I posted above with 01m data. Those two examples, one assigning CRS using gdal_edit.py, the other using gmt grdedit. I’ve got exactly the same results, so not posting anything here. gdal_edit.py assigned CRS correctly, gmt grdedit quietly failed

gmt grdimage @earth_gebco_01m_p -n+a -Aworld.tiff -J4326 produces error messages like below which seems like gdal errors (gdal seems to be called by grdimage internally to produce the .tiff?). However, the resulting images are totally fine.

ERROR 10: Pointer 'hTransform' is NULL in 'OCTTransform'.

ERROR 10: Pointer 'hTransform' is NULL in 'OCTTransform'.

ERROR 10: Pointer 'hTransform' is NULL in 'OCTTransform'.

More than 1000 errors or warnings have been reported. No more will be reported from now.

Yes, I can confirm the troubles when using tiled grids. This silently fails for me, which normally means a crash

gmt grdimage @earth_relief_05m -JX14/0 -P > lixo.ps       # Fine
psconvert lixo.ps -W+g        # Silently fails

Most weird because the PS file is ok. Note that here I used the GMT syntax -JX14/0 because I found another problem. While -J4326 works but assigns no CRS, using a proj4 string does. i.,e, -J+proj=lonlat. But in both cases it produces a square image and not even -J+proj=lonlat'width=14/0 changes that fact.

Current situation indicates that grdimage -A + gdal_edit is your best bet.

Found why grdimage -W+g fails for tilled grids. The PostScript lacks the necessary +proj=... information. Looks like when using tiles something is failing to add that info.

Note that -JEPSG or -J+proj are normally not supported in grdimage that most of times creates frames too (that’s where proj fails) but without -B that is not the case and it should work … but doesn’t.

However, the GMT syntax works. That is, this works

gmt grdimage @earth_relief_05m -JX14d/0d -P > lixo.ps
gmt psconvert lixo.ps -W+g

but still lacks the CRS info.

Thanks, it has been working for me at lower resolutions, but when I try to run it with a higher resolution image/dataset, the tile download fails part way through the set of tiles. I’m not sure why. I figure I can try to download the entire grd file from the server & use it as a local data source. This might also fix the error I get with grdimage & tiled downloads saving directly as TIFF.

grdblend [NOTICE]: Remote data courtesy of GMT data server oceania [http://oceania.generic-mapping-tools.org]

grdblend [NOTICE]: Earth Relief at 3x3 arc seconds tiles provided by SRTMGL3 (land only) [NASA/USGS].
grdblend [NOTICE]: → Download 1x1 degree grid tile (earth_relief_03s_g): S28E142
grdblend [ERROR]: Libcurl Error: Couldn’t resolve host name
grdblend [ERROR]: Unable to obtain remote file @S28E142.earth_relief_03s_g.nc
grdblend [NOTICE]: → Download 1x1 degree grid tile (earth_relief_03s_g): S28E142
grdblend [ERROR]: Libcurl Error: Couldn’t resolve host name
grdblend [ERROR]: Unable to obtain remote file @S28E142.earth_relief_03s_g.nc
grdblend [NOTICE]: → Download 1x1 degree grid tile (earth_relief_03s_g): S28E142
grdblend [ERROR]: Libcurl Error: Couldn’t resolve host name
grdblend [ERROR]: Unable to obtain remote file @S28E142.earth_relief_03s_g.nc
grdblend [ERROR]: File @S28E142.earth_relief_03s_g.nc not found
[Session gmt (0)]: Error returned from GMT API: GMT_FILE_NOT_FOUND (16)
grdimage [ERROR]: ERROR - Unable to produce blended grid from /tmp/=tiled_126_PO_cmYEAC
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
psconvert [ERROR]: world.ps: GMT PS format detected but file is not finalized. Maybe a -K in excess? No output created.

Yes, pre-downloading the tiles and assembling them locally will solve the current psconvert -W+g issue. I have a fix for it but needs more work.

Regarding the libcurl errors, don’t know, but you can try a mirror. See GMT_DATA_SERVER and Mirrors — The Generic Mapping Tools

Hi Joaqim,

I have rsynch’ed the entire server data content locally, but I can’t find info on how the tiled jp2 files are assembled into a whole grd file or image, can you point me at where some instructions are available?

Thanks…

Ah, that’s not how one do it. See gmtwhich -Gc. It will convert the the .jp2 into netCDF and store them in the cache dir (~./gmt/cache). grdcut -R -G... downloads only the necessary tiles (or use the ones in cache) and assembles a grid. Don’t know right now how to use the .jp2 that you already downloaded. Maybe put them in a html server and making GMT_DATA_SERVER point to it.

See Address issues raised in #8438 and more. by joa-quim · Pull Request #8440 · GenericMappingTools/gmt · GitHub