Something else must be going on here. Everything I tried with those .tiff shows correct locations and never a y with row numbers.
FYI, First cut at psconvert has a similar issue - but the lats seem sensible & the longs are way offset… Will persevere…
So something is screwed in my setup? Sigh… That makes it tricky if you can’t replicate it… especially with the brand new install of everything here.
Not sure what my next step is - it was all supposed to be relatively easy!!!
assigning CRS using gdal_edit.py works for me:
gmt grdimage @earth_gebco_10m_p -n+a -Aworld.tiff -J4326
gdal_edit.py -a_srs EPSG:4326 world.tiff
gmt begin quick png
gmt grdimage world.tiff
gmt coast -W -BSWen -Bxa -Bya --MAP_FRAME_TYPE=plain
gmt end show
using gmt grdedit quietly fails to assign CRS expressed as -J4326 to world.tiff, but the coastlines appear where they should:
gmt grdimage @earth_gebco_10m_p -n+a -Aworld.tiff -J4326
gmt grdedit world.tiff -J4326
gmt begin quick png
gmt grdimage world.tiff
gmt coast -W -BSWen -Bxa -Bya --MAP_FRAME_TYPE=plain
gmt end show
PS my gmt version is 6.5.0, part of my pygmt conda/mamba environment.
To be confirmed but I think grdedit works only with netCDF files.
Don’t know. Maybe install GMT6.5?
Can’t you reproduce this?
C:\v>gmt grdimage @earth_gebco_10m_p -J4326 -P > lixo.ps
C:\v>gmt psconvert lixo.ps -W+g
Input file size is 1654, 1654
0...10...20...30...40...50...60...70...80...90...100 - done.
C:\v>gdalinfo lixo.tiff
Driver: GTiff/GeoTIFF
Files: lixo.tiff
Size is 1654, 1654
Coordinate System is:
GEOGCRS["unknown",
DATUM["unknown",
ELLIPSOID["unknown",6378137,298.257220143428,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.217654171704510,90.000000000000000)
Pixel Size = (0.217654171705000,-0.108827085852000)
Metadata:
AREA_OR_POINT=Area
Software=GPL Ghostscript GIT PRERELEASE 9.54.0
TIFFTAG_XRESOLUTION=300
TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-180.2176542, 90.0000000) (180d13' 3.56"W, 90d 0' 0.00"N)
Lower Left (-180.2176542, -90.0000000) (180d13' 3.56"W, 90d 0' 0.00"S)
Upper Right ( 179.7823458, 90.0000000) (179d46'56.44"E, 90d 0' 0.00"N)
Lower Right ( 179.7823458, -90.0000000) (179d46'56.44"E, 90d 0' 0.00"S)
Center ( -0.2176542, 0.0000000) ( 0d13' 3.56"W, 0d 0' 0.00"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
Firstly, thanks heaps for your help. It really is appreciated!!!
Will try what you suggest - however - I just built a Win10 VM on the Mint box. Installed Miniconda, GMT & GDAL.
Then generated a world.tiff exactly as on Mint Linux (without the GMT prefix). Copied it to the Mint system & opened in QGIS - it overlays perfectly with GSHHG coastline.
So something is screwy with the Mint/UbuntuGIS setup… not really a GMT issue if it works perfectly on Windows.
Windows is not my preferred platform, but it seems to work. At least I can now do something useful!!!
I’m not sure if installing Miniconda on Linux, then gdal in Miniconda & running GMT from the Miniconda shell instead of a bash shell would generate correct geotiffs…
Thanks again!!!
And I just tried your grdimage/psconvert as above - that also worked, on Linux natively. Gave an image that overlays with QGIS.
If I then run: gdal_edit.py -a_crs EPSG:4326 test.tiff
The geotiff then has the correct CRS metadata embedded, as well as the extents.
Which sets me up (I think) to write a script to create the images/tiles I need…
Thanks again!!!
works for me
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.