That would be an awesome addition indeed 
From the gdal command given in the user guide, gdal_translate -of EHdr [..], it seems to be the EHdr driver (https://gdal.org/drivers/raster/ehdr.html). The AIG driver in gdal unfortunately seems not to support creation of datasets (see https://gdal.org/drivers/raster/index.html).
I’m trying this, but getting bad grids:
gmt grdcut -Gno.nc -R0/10/50/60 @earth_relief_05m
gdal_translate -of EHdr -ot Float32 no.nc no.flt
gmt grdinfo no.flt
no.flt: Title: no.hdr
no.flt: Command:
no.flt: Remark:
no.flt: Gridline node registration used [Cartesian grid]
no.flt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.flt: x_min: NaN x_max: NaN x_inc: NaN name: x n_columns: 0
no.flt: y_min: NaN y_max: NaN y_inc: NaN name: y n_rows: 0
no.flt: v_min: NaN v_max: NaN name: z
no.flt: scale_factor: 1 add_offset: 0
Done via gmt:
gmt grdcut -Gno.flt=gd:EHdr -R0/10/50/60 @earth_relief_05m
gmt grdinfo no.flt
no.flt: Title: no.hdr
no.flt: Command:
no.flt: Remark:
no.flt: Gridline node registration used [Cartesian grid]
no.flt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.flt: x_min: NaN x_max: NaN x_inc: NaN name: x n_columns: 0
no.flt: y_min: NaN y_max: NaN y_inc: NaN name: y n_rows: 0
no.flt: v_min: NaN v_max: NaN name: z
no.flt: scale_factor: 1 add_offset: 0
Guess it’s a gdal issue?
I’m getting side tracked now but;
Note format:
gmt grdcut -Gno.asc=ei -R0/10/50/60 @earth_relief_05m
grdcut [WARNING]: ESRI Arc/Info ASCII Interchange file must use proxy for NaN; default to -9999
gmt grdinfo no.asc
no.asc: Title:
no.asc: Command:
no.asc: Remark:
no.asc: Pixel node registration used [Cartesian grid]
no.asc: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.asc: x_min: 0 x_max: 10 x_inc: 0.0833333333333 name: x n_columns: 120
no.asc: y_min: 50 y_max: 60 y_inc: 0.0833333333333 name: y n_rows: 120
no.asc: v_min: NaN v_max: NaN name: z
no.asc: scale_factor: 1 add_offset: 0
v_min and v_max is both reported as NaN, but if you cat no.asc it contains only non-NaN’s.
Worked fine for me
grdcut @srtm_relief_03s -Ghawaii_south.flt=gd:EHdr -R-156/-155.5/19/19.5
gdalinfo hawaii_south.flt
Driver: EHdr/ESRI .hdr Labelled
Files: hawaii_south.flt
hawaii_south.flt.aux.xml
hawaii_south.hdr
hawaii_south.prj
Size is 601, 601
Coordinate System is:
...
But on a second look, it clipped the negative depths and replaced them by 0!!!
Origin = (-156.000416666666666,19.500416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left (-156.0004167, 19.5004167) (156d 0' 1.50"W, 19d30' 1.50"N)
Lower Left (-156.0004167, 18.9995833) (156d 0' 1.50"W, 18d59'58.50"N)
Upper Right (-155.4995833, 19.5004167) (155d29'58.50"W, 19d30' 1.50"N)
Lower Right (-155.4995833, 18.9995833) (155d29'58.50"W, 18d59'58.50"N)
Center (-155.7500000, 19.2500000) (155d45' 0.00"W, 19d15' 0.00"N)
Band 1 Block=601x1 Type=Float32, ColorInterp=Undefined
Min=-1.000 Max=4160.000
Minimum=-1.000, Maximum=4160.000, Mean=1172.830, StdDev=1136.810
NoData Value=nan
Metadata:
STATISTICS_MAXIMUM=4160
STATISTICS_MEAN=1172.8304905025
STATISTICS_MINIMUM=-1
STATISTICS_STDDEV=1136.809938332
STATISTICS_VALID_PERCENT=100
EDIT: Nope, it’s correct. The remote a SRTM has no negative nodes. We get the same if grdcut writes a nc file.
Ok, so the file is correct, but gmt is not able to report the ranges(?)
gmt grdcut -Gno.flt=gd:EHdr -R0/10/50/60 @earth_relief_05m
gmt grdinfo no.flt
no.flt: Title: no.hdr
no.flt: Command:
no.flt: Remark:
no.flt: Gridline node registration used [Cartesian grid]
no.flt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.flt: x_min: NaN x_max: NaN x_inc: NaN name: x n_columns: 0
no.flt: y_min: NaN y_max: NaN y_inc: NaN name: y n_rows: 0
no.flt: v_min: NaN v_max: NaN name: z
no.flt: scale_factor: 1 add_offset: 0
gdalinfo no.flt
Driver: EHdr/ESRI .hdr Labelled
Files: no.flt
no.flt.aux.xml
no.hdr
no.prj
Size is 120, 120
Coordinate System is:
GEOGCRS["GCS_unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["Degree",0.0174532925199433]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["Degree",0.0174532925199433]]]
Data axis to CRS axis mapping: 1,2
Origin = (0.000000000000000,59.999999999999964)
Pixel Size = (0.083333333333333,-0.083333333333333)
Corner Coordinates:
Upper Left ( 0.0000000, 60.0000000) ( 0d 0' 0.00"E, 60d 0' 0.00"N)
Lower Left ( 0.0000000, 50.0000000) ( 0d 0' 0.00"E, 50d 0' 0.00"N)
Upper Right ( 10.0000000, 60.0000000) ( 10d 0' 0.00"E, 60d 0' 0.00"N)
Lower Right ( 10.0000000, 50.0000000) ( 10d 0' 0.00"E, 50d 0' 0.00"N)
Center ( 5.0000000, 55.0000000) ( 5d 0' 0.00"E, 55d 0' 0.00"N)
Band 1 Block=120x1 Type=Float32, ColorInterp=Undefined
Min=-676.500 Max=1497.000
Minimum=-676.500, Maximum=1497.000, Mean=20.463, StdDev=217.651
NoData Value=nan
Metadata:
STATISTICS_MAXIMUM=1497
STATISTICS_MEAN=20.463229166667
STATISTICS_MINIMUM=-676.5
STATISTICS_STDDEV=217.65120068201
STATISTICS_VALID_PERCENT=100
It seems GMT needs a little help figuring out that it must call gdal
grdinfo hawaii_south.flt=gd
hawaii_south.flt: Title: Grid imported via GDAL
hawaii_south.flt: Command:
hawaii_south.flt: Remark:
hawaii_south.flt: Gridline node registration used [Geographic grid]
hawaii_south.flt: Grid file format: gd = Import/export through GDAL
hawaii_south.flt: x_min: -156 x_max: -155.5 x_inc: 0.000833333333333 (3 sec) name: x n_columns: 601
hawaii_south.flt: y_min: 19 y_max: 19.5 y_inc: 0.000833333333333 (3 sec) name: y n_rows: 601
hawaii_south.flt: v_min: -1 v_max: 4160 name: z
hawaii_south.flt: scale_factor: 1 add_offset: 0
+proj=longlat +datum=WGS84 +no_defs
Ah, it picked a wrong driver
I can confirm – adding =gd does the trick.
As some of you may have already seen, I added a PR to include the Leland’s code. If you build that GMT branch then you can explore it in Julia with
G = gmtread("hawaii_south.grd");
G2 = GMT.texture_img(G);
hill = gdaldem(G, "hillshade", zfactor=4);
G3 = blendimg!(G2, hill, new=true);
imshow(G3, fmt=:png)
The blendimg!() function let us blend a color image and a uint8 shade (as computed by gdaldem) using the equation shown above by Andreas or two uint8 images with transparency.
In Julia REPL just type ? blendimg! and ? GMT.texture_img(G) to see the online help.
Interested people please explore and create nice images that convince us that this is a good feature to add to plain GMT.
I’m waiting to see it too ! I’m interested to try in plain GMT
You can find @Joaquim’s work over at GitHub.
… which may or may not happen 
Some comparison between Traditional Hillshading and a combination with Texture Hillshading made with @Joaquim’s implementation – no colors:
The traditional way:
A combination of traditional hillshading and texture hillshading. Note the increased detail – especially in the shadows:
After staring too long at this stuff already – where do you see room for improvement @PlanetGus?
It’s already very impressive. Some colouring might help to highlight odd behaviours (on the alluvial fans for example ?)
How does the algorithm performs in this region ? -R10.13/10.43/45.46/45.85 (flat lands and high mountains)
Here is a quick test with your -R:
As you talked about alluvial fans I throw in a plot of the Death Valley National Park, CA, USA into the mix as well (Mesquite Flats area, -R-117.33/-116.84/36.51/36.84). Steep mountains and large alluvial fans:
And for those wanting color
G = gmtread("alps.grd", layout="TRB");
G2 = GMT.texture_img(G);
hill = gdaldem("alps.grd", "hillshade", zfactor=4);
color = gdaldem("alps.grd", "color-relief", color="faa.cpt");
G3 = blendimg!(G2, hill, new=true);
G4 = blendimg!(color, G3, new=true);
imshow(G4, fmt=:png)
This algorithm really gives a better visual perception. I like the bumps due to the alluvial fans 
Time to dig this up again after favourably talking about it in Oslo: Should I put this up on Github as a feature request?
Sorry. Yes, you can add a feature request with link to this page. Obviously not 6.5
Sorry, a feature request for what?
What I want is that branch merged.
1 Like
Sorry, Apri 21 is a while ago and I did not read this whole thread. Let me see. Have you merged master into that branch so it is ready?