How to generate a lat/long geotiff with grdimage

Brent · April 9, 2024, 12:42am

Hi,

I\m trying to generate an unprojected geotiff image with grdimage.

The image is a global relief map, including GEBCO & terrain data.

I get a good looking image, but can’t get the georeferencing correct. I want the coords to be in epsg:4326.

I’m trying things like:

gmt grdimage @earth_gebco_01m_p -Rd -Aworld.tiff -JEPSG:4326
gmt grdimage @earth_gebco_01m_p -R-180/180/-85/85 -Aworld.tiff -J"+proj=lonlat +datum=WGS84" +no_defs

But the output image has no CRS and the Y coords seem to be pixel counts, not latitudes.

I have set GDAL_DATA to /usr/share/gdal/ which is where all the GDAL csv files seem to be located.

I’m getting the error:
ERROR 10: Pointer ‘hTransform’ is NULL in ‘OCTTransform’.
repeatedly, but I don’t know why. It seems to relate to a bad EPSG code…

Any suggestions appreciated…

Joaquim · April 9, 2024, 10:40am

No idea. It works for me. Did not try the @earth_gebco_01m_p because don’t want to download it but tried with @earth_relief_01m_p

# On Win
grdimage @earth_gebco_02m_p -Rd -Aworld.tiff -JEPSG:4326
grdimage @earth_relief_01m_p -Rd -Aworld.tiff -JEPSG:4326

# On WSL
gmt grdimage @earth_gebco_10m_p -Rd -Aworld.tiff -JEPSG:4326

An old GMT version?

EDIT: Yes, I don’t see a projection info when a do gdalinfo but can open the file files in Mirone and QGis

mkononets · April 9, 2024, 11:17am

as a possible workaround, gdal_translate and gdal_edit have an option to assign a crs to a grid, -a_srs

gmt grdedit can do it too, using -J... as far as I can understand.

Joaquim · April 9, 2024, 12:42pm

The whole situation is confuse and clearly need improvements. For example this (Mollweide) works and stores the SRS, but it requires setting the map width when that is not needed.

gmt grdimage @earth_relief_10m_p -Aworld.tiff -JW15

And If I try with a proj4 string it errors. For map making we need to be able to convert from proj/WKT syntax in order to be able to create the frames, but with -A that is not the case so one should be able to pass the proj4 string or EPSG without doing a round trip to the GMT projection machinery.

Brent · April 9, 2024, 7:04pm

Thanks Joaquim,
I can open the output image in QGIS fine, but it is not correctly geo-located. Opening the GSHHG world coastline shapefile at the same time gives two maps on the canvas, they seem to be correctly aligned vertically, so the longitudes are somewhat consistent, but the latitudes are way different. The raster image is a small block at the top of the screen, the vector at the bottom.

I was hoping to use GMT to create a set of images/tiles for a global context layer served via WMS on our research vessels, but if the output images are not correctly georeferenced, it is not viable.

Brent

Brent · April 9, 2024, 7:06pm

Thanks for the suggestion, but it is not just the lack of a CRS, the Y coordinates are wrong as well, so just populating the CRS will not solve the problem.

Brent

Brent · April 9, 2024, 7:15pm

Thanks Joaquim,

Fresh install of GMT 6.4.4 on fresh Mint Linux setup on Ryzen based server- built from scratch today, then added UbuntuGIS unstable for the package installs via apt. GDAL is v3.8.4.

I had issues on my dated desktop so figured a completely new install from the OS up would be more robust…

I’m not sure why I keep getting the “ERROR 10: Pointer ‘hTransform’ is NULL…” message (or the first 1000 messages at least). Googling tells me it is a faulty EPSG code, but that does not seem to be the case.

Is this a possible clue to another issue on my setup?

Thanks,

Brent

Joaquim · April 9, 2024, 7:32pm

Odd, this is what I get with the Julia wrapper (pure GMT execution, just under other clothes).
Some stray warnings that don’t affect the result. Coastlines are plotted right.

julia> I = grdimage("@earth_relief_10m_p", A=true, Vd=1, B=:none, J=4326)
Warning: the following options were not consumed in grdimage => [:this_cpt]
        grdimage @earth_relief_10m_p  -J4326 -n+a -A
┌ Warning: Only 'I' for Images.jl and 'BRP' MEM layouts are allowed.
└ @ GMT C:\Users\j\.julia\dev\GMT\src\gmt_main.jl:486
A GMTimage object with 3 bands of type UInt8
Pixel node registration used
x_min: -180.0   x_max :180.0    x_inc :0.16666666666666666      n_columns :2160
y_min: -90.0    y_max :90.0     y_inc :0.16666666666666666      n_rows :1080
z_min: NaN      z_max :NaN
Mem layout:     TRPa

julia> viz(I, coast=true)
GMT [WARNING]: gmtapi_change_imagelayout: reordering function for case TRPa -> BRP not yet written. Doing nothing.

Joaquim · April 9, 2024, 7:34pm

And you get that error when doing what?

Remember that you can also create GeoTIFF images with psconvert -W

Brent · April 9, 2024, 7:52pm

I just ran essentially the same command:
gmt grdimage @earth_relief_10m_p -J4326 -n+a -Anext.tiff

This is the QGIS canvas… the cursor positions are fine for the shapefile, the image is offset.

Brent · April 9, 2024, 7:53pm

Yep, I will try psconvert as another approach…

Thanks!!

Joaquim · April 9, 2024, 8:09pm

Something else must be going on here. Everything I tried with those .tiff shows correct locations and never a y with row numbers.

Brent · April 9, 2024, 8:14pm

FYI, First cut at psconvert has a similar issue - but the lats seem sensible & the longs are way offset… Will persevere…

Brent · April 9, 2024, 8:17pm

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!!!

mkononets · April 9, 2024, 9:13pm

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.

Joaquim · April 9, 2024, 9:32pm

To be confirmed but I think grdedit works only with netCDF files.

Joaquim · April 9, 2024, 9:46pm

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

Brent · April 9, 2024, 10:35pm

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!!!

Brent · April 9, 2024, 10:47pm

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!!!

mkononets · April 9, 2024, 11:28pm

works for me