I’m a Julia beginner and try to figure out how to get imshow to color my terrain:
imshow(G3, fmt=:png, color="faa.cpt")
Doesn’t work for me I end up with a grayscale image and the following message:
Warning: the following options were not consumed in grdimage => [:color]
How do I correctly specify a CPT?
Hmm, it should have worked. Must investigate why it didn’t.
Notice, imshow is wrapper to grdimage and grdview and it does lots of guessing. Apparently some wizard power is failing in this case. I’m betting that’s because you are trying to plot an image and setting the cpt, which probably GMT itself does not accept. Now that I’m writing this it does ring a bell.
As a general advise to learn the GMT,jl way, use the option Vd to print the GMT command that is generated under the hood. Vd=1 prints and executes, Vd=2 just prints it. E.g.
imshow(G3, color="faa.cpt", Vd=2)
Warning: the following options were not consumed in grdimage => [:color]
"grdimage -R10.5098611111111/11.4101388888889/46.6298611111111/47.1401388888889 -JX14c/0 -Baf -BWSen -n+a -D -P -K > C:\\TMP\\GMTjl_tmp.ps"
The other advise is to read at least this
Yes, that’s on purpose. You cannot set a CPT in grdimage when input is an image.
But this case should caught earlier and print a meaningful error message. Please open an issue in GMT.jl for this
Besides, the solution for this case is failing too. And that would be to use the image_cpt! function. This should have worked
C = gmtread("faa.cpt");
image_cpt!(G3, C)
imshow(G3)
Wise words … I skimmed it but more in-depth reading & learning is definitely necessary. Thank you for the Vd=1 and Vd=2 tips!
Is that the error you are talking about?
G = gmtread("alps2.grd", layout="TRB");
G2 = GMT.texture_img(G);
makecpt(cmap=:dem4, cptname="my_dem4.cpt");
C = gmtread("my_dem4.cpt");
hill = gdaldem("alps2.grd", "hillshade", zfactor=2.5);
color = gdaldem("alps2.grd", "color-relief", C);
G3 = blendimg!(G2, hill, new=true, transparency=0.61);
G4 = blendimg!(color, G3, new=true);
imshow(G4, fmt=:png)
ERROR: LoadError: MethodError: no method matching helper_run_GDAL_fun(::typeof(gdaldem), ::String, ::String, ::GMT.GMTcpt, ::String)
Nope. The error in my example is much more illusive. Everything seems to run fine but the produced PS is corrupted. I’m afraid I’ll have to dig into the GMT guts to see what is failing.
Note that these two
makecpt(cmap=:dem4, cptname="my_dem4.cpt");
C = gmtread("my_dem4.cpt");
can be collapsed in
C = makecpt(cmap=:dem4);
All modules that produce files can be used to save that file in disk or just return the data to be used in Julia, wrapped in the appropriate type
Will look why your example fail.
Ok, this is the current version:
using GMT
Gt = gmtread("alps2.grd", layout="TRB");
Gh = gmtread("alps2.grd");
G2 = GMT.texture_img(Gt);
C = makecpt(cmap=:dem4);
hill = gdaldem(Gh, "hillshade", zfactor=2.5);
color = gdaldem(Gh, "color-relief", C, Vd=1);
G3 = blendimg!(G2, hill, new=true, transparency=0.61);
G4 = blendimg!(color, G3, new=true);
imshow(G4, fmt=:png)
And this the full error message:
> julia alps.jl
ERROR: LoadError: MethodError: no method matching helper_run_GDAL_fun(::typeof(gdaldem), ::GMT.GMTgrid{Float32, 2}, ::String, ::GMT.GMTcpt, ::String, ::Pair{Symbol, Int64})
Closest candidates are:
helper_run_GDAL_fun(::Function, ::Any, ::String, ::Vector{String}, ::String, ::Any...) at /home/kristof/.julia/packages/GMT/dCcJn/src/gdal_tools.jl:90
helper_run_GDAL_fun(::Function, ::Any, ::String, ::Vector{String}) at /home/kristof/.julia/packages/GMT/dCcJn/src/gdal_tools.jl:90
Stacktrace:
[1] gdaldem(indata::GMT.GMTgrid{Float32, 2}, method::String, opts::GMT.GMTcpt; dest::String, kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:Vd,), Tuple{Int64}}})
@ GMT.Gdal ~/.julia/packages/GMT/dCcJn/src/gdal_tools.jl:85
[2] top-level scope
@ ~/Dokumente/GMT/teximg/alps.jl:7
in expression starting at /home/kristof/Dokumente/GMT/teximg/alps.jl:7
>
The line 7 it is complaining about is this one:
color = gdaldem(Gh, "color-relief", C, Vd=1);First, in my defense, I’m trying to make the use of gdaldem more user friendly. Specially in light of GDAL claims that they support GMT CPT’s. which is basically false.
I think I know where problem is. gdaldem(..., "hillshade", ...) was tested and use smoothed out but gdaldem(..., "color-relief", ...) much less. I think that if you use color="cptfilename.cpt" it will work but passing the GMTcpt structure is not being recognized (or some shit around that)
Last, The Vd mechanism works only for the GMT modules calls, not for GDAL’s (different machinery and solution, if possible, not yet invented).
No need for defensive actions – I tip my hat for your work that went into GMT.jl and your ongoing support. Thank you! 
color="cptfilename.cpt" does indeed work.
Not that happy with the blending of the colored terrain and the hillshading and for some reason I can’t use blendimg!(color, G3, new=true, transparency=0.4); It seems to dislike the transparency part and throws an error.
No, you can’t use transparency when blending a color image and texture. It was not in the equation that Andreas posted. Nothing than cannot be changed but so far it’s like that.
Update your GMT.jl and you can now do
color = gdaldem("alps.grd", "color-relief", color=C);
To update the Julia wrapper you need to be using its master version. Don’t know if you have to uninstall the version you have (] rm GMT) before installing the master with
] add GMT#master
You may be able to mix in some transparency via grdmix.
That’s a great idea – thank you @pwessel. It is not really transparency what I’m after, it is more of a weighting/mixing thing.
Prior to my first Julia steps I did quite some modifications to the intensity grid file with grdmath like compressing the light areas (>0) into the (0, 0.1) range and the dark areas (<0) into the (-0.7, 0) range.
I’m hoping that you and the core devs see texture shading as a useful permanent addition to GMT leading to an easier integration into my existing scripts. That Julia language seems to be a quite useful rabbit hole to explore …
The hole texture thing (excluding the blending) in Julia boils down to
function texture_img(G::GMTgrid; detail=1.0, contrast=2.0, uint16=false)
texture = deepcopy(G.z)
terrain_filter(texture, detail, size(G,1), size(G,2), G.inc[1], G.inc[2], 0)
(startswith(G.proj4, "+proj=merc")) && fix_mercator(texture, detail, size(G,1), size(G,2), G.range[3], G.range[4])
terrain_image_data(texture, contrast, size(G,1), size(G,2), 0.0, (uint16) ? 65535.0 : 255.0)
mat = (uint16) ? round.(UInt16, texture) : round.(UInt8, texture)
Go = mat2img(mat, hdr=grid2pix(G), proj4=G.proj4, wkt=G.wkt, noconv=true, layout=G.layout*"a")
Go.range[5:6] .= extrema(Go.image)
Go
end
with these two wrapper functions to talk directly with the lib.
This is why I think Julia brings another dimension to GMT expansion
terrain_filter(data, detail, nrows, ncols, xinc, yinc, coord_type, center_lat=0.0, progress=C_NULL) =
ccall((:terrain_filter, thelib), Cint, (Ptr{Cfloat}, Cdouble, Cint, Cint, Cdouble, Cdouble, Cint, Cdouble, Ptr{Cvoid}), data, detail, nrows, ncols, xinc, yinc, coord_type, center_lat, progress)
terrain_image_data(data, contrast, nrows, ncols, image_min=0., image_max=65535.0) =
ccall((:terrain_image_data, thelib), Cint, (Ptr{Cfloat}, Cint, Cint, Cdouble, Cdouble, Cdouble), data, nrows, ncols, contrast, image_min, image_max)