From 9c8a4dd5f5445c1354911d212086543167047a84 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 10 Feb 2024 18:43:45 +0000 Subject: [PATCH 1/5] slicecube() accepts subregions. --- src/utils_types.jl | 81 +++++++++++++++++++++++++++++----------------- 1 file changed, 52 insertions(+), 29 deletions(-) diff --git a/src/utils_types.jl b/src/utils_types.jl index efdb6c7fd..cdc59ebe1 100644 --- a/src/utils_types.jl +++ b/src/utils_types.jl @@ -1196,7 +1196,7 @@ Extract a slice from a GMTgrid cube. neighboring layers. For example, `slice=2.5` on a cube were layers are one unit apart will interpolate between layers 2 and 3 where each layer weights 50% in the end result. NOTE: the return type is still a cube but with one layer only (and the corresponding axis coordinate). - slice` Can also be a vector of integers representing the slices we want to extract. The output is another cube. + `slice` Can also be a vector of integers representing the slices we want to extract. The output is another cube. - `axis`: denotes the dimension being sliced. The default, "z", means the slices are taken from the vertical axis. `axis="x"` means slice along a column, and `axis="y"` slice along a row. @@ -1238,8 +1238,10 @@ function slicecube(I::GMTimage, layer::Union{Int, AbstractVector{<:Int}}) GMTimage(I.proj4, I.wkt, I.epsg, I.geog, range, copy(I.inc), I.registration, I.nodata, "Gray", I.metadata, names, copy(I.x), copy(I.y), [0.], mat, zeros(Int32,3), String[], 0, Array{UInt8,2}(undef,1,1), I.layout, I.pad) end -function slicecube(G::GMTgrid, slice::Union{Int, AbstractVector{<:Int}}; x::VecOrMat{<:Real}=Float64[], y::VecOrMat{<:Real}=Float64[], axis="z") +function slicecube(G::GMTgrid, slice::Union{Int, AbstractVector{<:Int}}; axis="z", + x::Union{VecOrMat{<:Real}, Tuple}=Float64[], y::Union{VecOrMat{<:Real}, Tuple}=Float64[]) # Method that slices grid cubes. SLICE is the row|col|layer number. AXIS picks the axis to be sliced + # This method lets us slice a cube along any or all of the x|y|z axis (ndims(G) < 3 || size(G,3) < 2) && error("This is not a cube grid.") _axis = lowercase(string(axis)) @@ -1252,33 +1254,48 @@ function slicecube(G::GMTgrid, slice::Union{Int, AbstractVector{<:Int}}; x::VecO isempty(G.v) && (G.v = collect(1:size(G,3))) rng = isvec ? slice : slice:slice - colmajor = (G.layout[2] == 'C') - if (_axis == "z") - G_ = mat2grid(G[:,:,slice], G.x, G.y, isvec ? G.v[slice] : [G.v[slice]], reg=G.registration, is_transposed=!colmajor) + colmajor, ix, iy, ixp, iyp = helper_slicecube(G, x, y) + + if (_axis == "z") # A horizontal slice (plane xy) + _ix, _iy = colmajor ? (ix, iy) : (iy, ix) + G_ = mat2grid(G[_iy,_ix,slice], G.x[ixp], G.y[iyp], isvec ? G.v[slice] : [G.v[slice]], reg=G.registration, is_transposed=!colmajor) G_.names = (!isempty(G.names) && !all(G.names .== "")) ? (isvec ? G.names[slice] : [G.names[slice]]) : G.names G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = G.proj4, G.wkt, G.epsg, G.geog, G.layout - elseif (_axis == "y") + elseif (_axis == "y") # A slice in xz plane + (!colmajor && G.layout[1] == 'T') && (rng = (size(G,2)+1 - rng[end]):(size(G,2)+1 - rng[1])) # Believe that (10 - rng) errors but (10 - 1:1) NOT!!!!! yv = isvec ? G.y[slice] : [G.y[slice]] - t = isempty(x) ? [1 colmajor ? size(G,2) : size(G,1)] : round.(Int,interp_vec(G.x, x)) - ix = t[1]:t[2] - if (colmajor) G_ = mat2grid(G[rng,ix,:], G.x[ix], yv, G.v, reg=G.registration, names=G.names) - else G_ = mat2grid(G[ix,rng,:], G.x[ix], yv, G.v, reg=G.registration, is_transposed=true, names=G.names) + if (colmajor) G_ = mat2grid(G[rng,ix,:], G.x[ixp], yv, G.v, reg=G.registration, names=G.names) + else G_ = mat2grid(G[ix,rng,:], G.x[ixp], yv, G.v, reg=G.registration, is_transposed=true, names=G.names) end G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = "", "", 0, 0, "TRB" - else + else # A slice in yz plane + _iy = (!colmajor && G.layout[1] == 'T') ? ((size(G,2)+1 - iy[end]):(size(G,2)+1 - iy[1])) : iy xv = isvec ? G.x[slice] : [G.x[slice]] - t = isempty(y) ? [1 colmajor ? size(G,1) : size(G,2)] : round.(Int,interp_vec(G.y, y)) - iy = t[1]:t[2] - if (colmajor) G_ = mat2grid(G[iy,rng,:], xv, G.y[iy], G.v, reg=G.registration, names=G.names) - else G_ = mat2grid(G[rng,iy,:], xv, G.y[iy], G.v, reg=G.registration, is_transposed=true, names=G.names) + if (colmajor) G_ = mat2grid(G[iy,rng,:], xv, G.y[iyp], G.v, reg=G.registration, names=G.names) + else G_ = mat2grid(G[rng,_iy,:], xv, G.y[iyp], G.v, reg=G.registration, is_transposed=true, names=G.names) + (G.layout[1] == 'T') && (G_.z = fliplr(G_.z)) # The debugger told me to do this. end G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = "", "", 0, 0, "TRB" end return G_ end -function slicecube(G::GMTgrid, slice::AbstractFloat; x::VecOrMat{<:Real}=Float64[], y::VecOrMat{<:Real}=Float64[], axis="z") +function helper_slicecube(G, x, y) + # Shared between two methods of slicecube. + colmajor = (G.layout[2] == 'C') + tx = isempty(x) ? [1 colmajor ? size(G,2) : size(G,1)] : round.(Int,interp_vec(G.x, x)) + ix = tx[1]:tx[2] + ty = isempty(y) ? [1 colmajor ? size(G,1) : size(G,2)] : round.(Int,interp_vec(G.y, y)) + iy = ty[1]:ty[2] + # The pixel registration case is more complicated. Not sure that the following is the right way to do it. + ixp, iyp = (G.registration == 0) ? (ix, iy) : (ix[1]:ix[end]+1, iy[1]:iy[end]+1) + return colmajor, ix, iy, ixp, iyp +end + +function slicecube(G::GMTgrid, slice::AbstractFloat; axis="z", + x::Union{VecOrMat{<:Real}, Tuple}=Float64[], y::Union{VecOrMat{<:Real}, Tuple}=Float64[]) # Method that slices grid cubes. SLICE is the x|y|z coordinate where to slice. AXIS picks the axis to be sliced + # So far horizontal slices are unique (single slice) but vertical slices can slice sub-cubes. (ndims(G) < 3 || size(G,3) < 2) && error("This is not a cube grid.") _axis = lowercase(string(axis)) @@ -1290,26 +1307,28 @@ function slicecube(G::GMTgrid, slice::AbstractFloat; x::VecOrMat{<:Real}=Float64 all(frac .< 0.1) && return slicecube(G, layer, axis=_axis, x=x, y=y) # If 'slice' is within 10% of lower or upper layer just all(frac .> 0.9) && return slicecube(G, layer.+1, axis=_axis, x=x, y=y) # return that layer and do not interpolation between layers. - nxy = size(G,1)*size(G,2) + nxy = size(G,1) * size(G,2) l = layer - colmajor = (G.layout[2] == 'C') - if (_axis == "z") + colmajor, ix, iy, ixp, iyp = helper_slicecube(G, x, y) + + if (_axis == "z") # A horizontal slice (plane xy) mat = [G[k] + (G[k+nxy] - G[k]) * frac for k = (layer-1)*nxy+1 : layer*nxy] G_ = mat2grid(reshape(mat,size(G,1),size(G,2)), G.x, G.y, [Float64(slice)], reg=G.registration, is_transposed=!colmajor, layout=G.layout) - elseif (_axis == "y") - t = isempty(x) ? [1 colmajor ? size(G,2) : size(G,1)] : round.(Int,interp_vec(G.x, x)) - ix = t[1]:t[2] + if (!isempty(x) || !isempty(y)) # It would have been too complicated to do this with the "mat = ..." above + isempty(x) && (x = [G.range[1], G.range[2]]) + isempty(y) && (y = [G.range[3], G.range[4]]) + G_ = crop(G_, region=(x[1], x[2], y[1], y[2])) + end + elseif (_axis == "y") # A slice in xz plane if (colmajor) mat = G[l:l,ix,:] .+ (G[l+1:l+1,ix,:] .- G[l:l,ix,:]) .* frac else mat = G[ix,l:l,:] .+ (G[ix,l+1:l+1,:] .- G[ix,l:l,:]) .* frac end - G_ = mat2grid(mat, G.x[ix], [G.y[l]], G.v, reg=G.registration, is_transposed=!colmajor, layout="TRB") - else - t = isempty(y) ? [1 colmajor ? size(G,1) : size(G,2)] : round.(Int,interp_vec(G.y, y)) - iy = t[1]:t[2] + G_ = mat2grid(mat, G.x[ixp], [G.y[l]], G.v, reg=G.registration, is_transposed=!colmajor, layout="TRB") + else # A slice in yz plane if (colmajor) mat = G[iy,l:l,:] .+ (G[iy,l+1:l+1,:] .- G[iy,l:l,:]) .* frac else mat = G[l:l,iy,:] .+ (G[l+1:l+1,iy,:] .- G[l:l,iy,:]) .* frac end - G_ = mat2grid(mat, [G.x[l]], G.y[iy], G.v, reg=G.registration, is_transposed=!colmajor, layout="TRB") + G_ = mat2grid(mat, [G.x[l]], G.y[iyp], G.v, reg=G.registration, is_transposed=!colmajor, layout="TRB") end return G_ end @@ -1365,12 +1384,16 @@ const cubeslice = slicecube # I'm incapable of remembering which one it is. # --------------------------------------------------------------------------------------------------- """ + Gs = squeeze(G::GMTgrid) -> GMTgrid + +Remove singleton dimension from a grid. So far only for vertical slices of 3D grids. """ function squeeze(G::GMTgrid) dims = size(G) if ((ind = findfirst(dims .== 1)) !== nothing) - (ind == 1) && return mat2grid(reshape(G.z, dims[2], dims[3]), x=G.x, y=G.v, layout="TRB", is_transposed=true) - (ind == 2) && return mat2grid(reshape(G.z, dims[1], dims[3]), x=G.y, y=G.v, layout="TRB", is_transposed=true) + which_x = (length(G.x) > 1) ? G.x : G.y # Rather than relying in mem_layout, pick the non-singleton vector + (ind == 1) && return mat2grid(reshape(G.z, dims[2], dims[3]), x=which_x, y=G.v, layout="TRB", is_transposed=true) + (ind == 2) && return mat2grid(reshape(G.z, dims[1], dims[3]), x=which_x, y=G.v, layout="TRB", is_transposed=true) end end From e819ee8ce37089bb02599fa86c0488348286bf48 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 10 Feb 2024 18:45:18 +0000 Subject: [PATCH 2/5] cubeplot() now accepts inset and received a new method to work with cubes. --- src/utils_project.jl | 222 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 212 insertions(+), 10 deletions(-) diff --git a/src/utils_project.jl b/src/utils_project.jl index f8ab33fac..c7fc0920f 100644 --- a/src/utils_project.jl +++ b/src/utils_project.jl @@ -343,10 +343,10 @@ Plot grid lines on top of an image created with the `worldrectangular` function. - `GI`: A GMTgrid or GMTimage data type. - `grid`: A vector of GMTdatset with meridians and parallels to be plotted. This is normaly produced - by the `graticules()` or `worldrectgrid()` functions. + by the `graticules()` or `worldrectgrid()` functions. - `annot`: Wether to plot coordinate annotations or not (`annot=false`). - `sides`: Which sides of plot to annotate. `W` or `L` means annotate the left side and so on for any - combination of "WESNLRBT". To not annotate a particular side just omit that character. *e.g.* + combination of "WESNLRBT". To not annotate a particular side just omit that character. *e.g.* `sides="WS"` will annotate only the left and bottom axes. - `figname`: To create a figure in local directory and with a name `figname`. If `figname` has an extension that is used to select the fig format. *e.g.* `figname=fig.pdf` creates a PDF file localy called 'fig.pdf' @@ -430,7 +430,7 @@ end # ----------------------------------------------------------------------------------------------- """ cubeplot(img1::Union{GMTimage, String}, img2::Union{GMTimage, String}="", img3::Union{GMTimage, String}=""; - back::Bool=false, show=false, notop::Bool=false, kw...) + back::Bool=false, show=false, notop::Bool=false, xlabel="", ylabel="", zlabel="", title="", kw...) Plot images on the sides of a cube. Those images can be provided as file names, or GMTimage objects. @@ -446,19 +446,25 @@ Plot images on the sides of a cube. Those images can be provided as file names, The `kw...` keyword/value options may be used to pass: -- `region`: The limits extents that will be used to annotate the *x,y,z* axes. It uses the same syntax as all +- `region, limits`: The limits extents that will be used to annotate the *x,y,z* axes. It uses the same syntax as all other modules that accept this option (*e.g.* ``coast``). It defaults to "0/9/0/9/-9/0" - `figsize`: Select the horizontal size(s). Defaults to 15x15 cm. - `zsize`: Sets the size of *z* axis in cm. The default is 15. - `view`: The view point. Default is `(135,30)`. WARNING: only azimute views from the 4rth quadrant are implemented. - `transparency`: Sets the image's transparency level in the [0,1] or [0 100] intervals. Default is opaque. - `inset` or `hole`: Draws an inset hole in the cube's Southern wall. This option is a tuple with the form: - ``((img1,img2[,img3]), width=?, [depth=?])`` where ``(img1,img2[,img3])`` are the images of the `west`, `north` + ``((img1,img2[,img3]), width=?, [depth=?])`` where ``(img1,img2[,img3])`` are the images of the `north` + (that is, the plane whose normal is `y`) and `west` (that is, the plane whose normal is `x`) and, optionally, `bottom` sides of the inset. The `width` and `depth` are the width and depth of the the inset. If `depth` is not provided it defauls to `width`. These values **must** be given in percentage of the cube's width and can be given in the [0-1] or [0-100] interval. +- `xlabel, ylabel, zlabel, title`: Optional axes labels and title. Each one of these must be a string. +- `cmap, colormap, cpt, colorscale`: Add a colorbar at the bottom of the figure. The colormap can be passed as a + single argument or as a tuple of arguments, where first must be the colormap and second [and optional a third] + are the axes colorbar labels taking the form: `cmap=(C, "xlabel=Blabla1"[, "ylabel=Blabla2"])`. """ -function cubeplot(fname1::Union{GMTimage, String}, fname2::Union{GMTimage, String}="", fname3::Union{GMTimage, String}=""; back::Bool=false, notop::Bool=false, show=false, kw...) +function cubeplot(fname1::Union{GMTimage, String}, fname2::Union{GMTimage, String}="", fname3::Union{GMTimage, String}=""; + back::Bool=false, notop::Bool=false, xlabel="", ylabel="", zlabel="", title="", show=false, kw...) # ... d = KW(kw) opt_R = ((txt::String = parse_R(d, "")[2]) != "") ? txt[4:end] : "0/9/0/9/-9/0" @@ -475,18 +481,32 @@ function cubeplot(fname1::Union{GMTimage, String}, fname2::Union{GMTimage, Strin else f2, f3 = fname2, fname2 # Two images, repeat second on the vert sides end - basemap(R=opt_R, J=opt_J, JZ=opt_JZ, p=opt_p) + basemap(R=opt_R, J=opt_J, JZ=opt_JZ, p=opt_p, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, title=title) + if ((val = find_in_dict(d, CPTaliases)[1]) !== nothing) # See if want to plot a colorbar. + if isa(val, GMTcpt) + colorbar!(val, position=(outside=true, anchor=:BC, offset="0.9c"), B=:af) + elseif isa(val, Tuple) + len::Int = length(val) + C::GMTcpt = val[1] # If first arg is not a CPT, an error occurs here. + xlabel = startswith(val[2], "xlabel=") ? val[2][8:end] : "" # These 4 lines are to allow both combinations + ylabel = startswith(val[2], "ylabel=") ? val[2][8:end] : "" + (xlabel == "" && len == 3) && (xlabel = startswith(val[3], "xlabel=") ? val[3][8:end] : "") + (ylabel == "" && len == 3) && (ylabel = startswith(val[3], "ylabel=") ? val[3][8:end] : "") + colorbar!(C, position=(outside=true, anchor=:BC, offset="0.9c"), B=:af, xlabel=xlabel, ylabel=ylabel) + end + end + vsz = parse(Float64, opt_JZ) TB = front ? :T : :B SN = front ? :S : :N EW = front ? :E : :W bak = CTRL.pocket_J[3] # Save this because sideplot() calls parse_JZ with too few info to preserve it in case of need. - !notop && image!(f1, compact=sideplot(plane=TB, vsize=vsz), t=opt_t) # 'compact' is that option that lets pass a GMT str + !notop && image!(f1, compact=sideplot(plane=TB, vsize=abs(vsz)), t=opt_t) # 'compact' is that option that lets pass a GMT str image!(f2, compact=sideplot(plane=SN, vsize=vsz), t=opt_t) R = image!(f3, compact=sideplot(plane=EW, vsize=vsz), t=opt_t) CTRL.pocket_J[3] = bak - if ((val = find_in_dict(d, [:hole :inset])[1]) !== nothing) + if ((val = find_in_dict(d, [:hole :inset])[1]) !== nothing) # If we have an inset (hole) request (!isa(val, Tuple) || (!isa(val[1], Tuple) && 2 <= length(val) <= 3)) && error("The 'inset' option must be a tuple with the form ((img1,img2[,img3]), width=?, [depth=?])") hole_width::Float64 = val[2] @@ -498,7 +518,9 @@ function cubeplot(fname1::Union{GMTimage, String}, fname2::Union{GMTimage, Strin azim = parse(Float64, split(opt_p, '/')[1]) width = parse(Float64, split(opt_J[2:end], '/')[1]) - hole_width, hole_depth = hole_width * width, hole_depth * width + spliR = parse.(Float64,split(opt_R, '/')) + height = width * (spliR[4] - spliR[3]) / (spliR[2] - spliR[1]) + hole_width, hole_depth = hole_width * width, hole_depth * height xoff = (width - hole_width) * cosd(azim -90) td = (hole_depth > 0.0) ? "$hole_depth" : "0" # Found that "/0.0" screws up basemap!(R=opt_R, J="X$(hole_width)/"*td, JZ=opt_JZ, p=opt_p, X=xoff, B=0) @@ -571,3 +593,183 @@ function sideplot(; plane=:xz, vsize=8, depth=NaN, kw...) _W, _H = (p == 'x') ? (H, zsize) : (p == 'y') ? (W, zsize) : (W, H) @sprintf(" -Dg%.12g/%.12g+w%.12g/%.12g -p%c%.0f/%.0f %s %s", lims[1], lims[3], _W, _H, p, azim, elev, opt_X, opt_Y) end + +# ----------------------------------------------------------------------------------------------- +""" + cubeplot(G::GMTgrid; top=nothing, topshade=false, zdown::Bool=false, xlabel="", ylabel="", zlabel="", title="", + show=false, interp::Float64=0.0, kw...) + +Make a 3D plot of a 3D GMTgrid (a cube) with a top view perspective from the 4rth quadrant (only one implemented). +There are several options to control the paintinf of the cube walls but off course not all possibilities are covered. +For ultimate control, users cab create the side wall images separately and feed them to the cubeplot method that +accepts only images as input. + +- `cmap, colormap, cpt, colorscale`: Pass in a GMTcpt colormap to be used to paint the vertical walls and + optionally, the top wall. The default is to compute this from the cube's min/max values with the `turbo` colormap. + +- `colorbar`: Add a colorbar at the bottom of the figure. The plotted colormap is either the auto-generated colormap + (from the cube's min/max and the `turbo` colormap) or the one passed via the `cmap` option. The optional syntax of + this option is either: `colorbar=true` or `colorbar=(C, "xlabel=Blabla1"[, "ylabel=Blabla2"])`. Attention that when + the labels request are passed, thy MUST conform with the `xlabel=...` and `ylabel=...` prefix part. + +- `inset`: Add an inset to the figure. This inset takes the form of a _hole_ located in the lower right corner of + the cube in which its inner walls are painted with partial vertical slices of the cube. The `inset` option + may be passed as a two elements array or tuple where first element is the satrting longitude (end is cube's eastermost + coordinate) and second the ending latitude (start is southernmost lat): an alternative syntax is to use + `inset=(lon=?, lat=?)`. + +- `interp`: When the cube layers are not equi-distant, the vertical side walls are not regular gris. This option, + that is called by default, takes care of obtaining a regular grid by linear interpolation along the columns. + This automatic interpolation uses the smallest increment in the vertical direction, but that may be overridden + by the `interp` increment option (a float value). + +- `region, limits`: A 4 elements array or Tuple (`x_min, x_max, y_min, y_max`) with the limits of a sub-region to display. + Default uses the entire cube. + +- `show`: If `true` display the figure. Default is `false`, _i.e._ it lets append further elements latter if wished. + +- `top`: An optional GMTgrid or the file name of a GMTgrid to be used to create the top wall of the cube. If, instead, + a GMTimage is passed we plot it directly (grids are converted to images using default colormaps or one passed via `topcmap`). + The default is to use the cube's first slice as the top wall. + +- `topshade`: Only used when the `top` option was used to pass a GMTgrid (or the name of one). When `true`, the top wall + image is created with the cube's first slice and a shaded effect computed with `grdgradient` on the `top` grid + (normally a topography grid). This creates a nice effect that shows both the cube's first layer and the topography where it lies. + Optionally, the `topshade` option may be used with a `grdgradient` grid computed with any other grid. + +- `topcmap`: An alternative colormap for the top wall. Default is the same as the `turbo` computed with cube `G` itself. + +- `xlabel, ylabel, zlabel, title`: Optional axes labels and title. Each one of these must be a string. + +- `zdown`: When `true`, the z-axis is positive down. Default is `false` (positive up). + +- `zsize`: Vertical size of the plotted cube. Default is 6 cm. Use a negative value if the z-axis is positive down, but + see also the alternative `zdown` option. + +### Example: +```julia +C = xyzw2cube("USTClitho2.0.wrst.sea_level.txt"); +cubeplot(C, top="@earth_relief", inset=(lon=110,y=40), topshade=true, zdown=true, title="Vp model", + colorbar=("xlabel=P-wave velocity","ylabel=km/s"), show=true) +``` +""" +function cubeplot(G::GMTgrid; top=nothing, topshade=false, zdown::Bool=false, xlabel="", ylabel="", zlabel="", title="", + show=false, interp::Float64=0.0, kw...) + (size(G,3) < 2) && error("First input must be a 3D grid.") + d = KW(kw) + if ((opt_R = parse_R(d, "", false)[2]) != "") + x_min, x_max, y_min, y_max, = opt_R2num(opt_R) + else + x_min, x_max, y_min, y_max = G.range[1], G.range[2], G.range[3], G.range[4] + end + zsize = ((opt_JZ = parse_JZ(d, "")[2]) != "") ? parse(Float64, opt_JZ[5:end]) : 6.0 + (zdown && zsize > 0) && (zsize *= -1) # This way we can use a negative size directly or use the zdown opt + ((C = find_in_dict(d, CPTaliases)[1]) === nothing) && (C = makecpt(G.range[5],G.range[6])) + ((Ctop = find_in_dict(d, [:topcmap])[1]) === nothing) && (Ctop = C) + showcpt, cbar_label = false, "" + if ((val = find_in_dict(d, [:colorbar])[1]) !== nothing) + showcpt = true + (isa(val, String) || isa(val, Tuple)) ? (cbar_label = val) : (!isa(val, Integer) && + error("The 'colorbar' option must be a string or a tuple with two strings.")) + end + + # See about the inset option + inset::Vector{Float64}=Float64[] + if ((val = find_in_dict(d, [:inset])[1]) !== nothing) + (isa(val, VecOrMat) || isa(val, Tuple)) && (inset = [val[1], val[2]]) + if (isa(val, NamedTuple)) + ks::Tuple{Symbol, Symbol} = keys(val) + (:lon in ks) ? append!(inset, val[:lon]) : ((:x in ks) ? append!(inset, val[:x]) : nothing) + (:lat in ks) ? append!(inset, val[:lat]) : ((:y in ks) ? append!(inset, val[:y]) : nothing) + (length(inset) != 2) && error("The 'inset' option must have to elements. The lon and lat limits of the inset.") + end + end + + # Decide what to put at the top of the pile + if (isa(top, String) || isa(top, GMTgrid)) + R = (x_min,x_max,y_min,y_max) + Gt = (isa(top, String)) ? ((top[1] == '@') ? grdcut(top, R=R, J="Q15") : gmtread(top, R=R)) : crop(top, R=R) + if (topshade == 1) + It = grdimage(grdsample(slicecube(G,1), R=Gt), A=true, B=:none, C=C, shade=grdgradient(Gt,A=-45,N=:t), layout="BRP") + else + It = grdimage(Gt, A=true, B=:none, C=:topo, shade=true, layout="BRP") + end + elseif (isa(top, GMTimage)) + It = crop(top, region=(x_min,x_max,y_min,y_max)) + else + It = grdimage(slicecube(G,1), A=true, B=:none, layout="BRP", C=C) + end + + Gs = interp_vslice(squeeze(slicecube(G, G.range[3], axis="y")), inc=interp) # The South wall + Is = grdimage(Gs, A=true, B=:none, layout="BRP", C=C) + Ge = interp_vslice(squeeze(slicecube(G, G.range[2], axis="x")), inc=interp) # The East wall + Ie = grdimage(Ge, A=true, B=:none, layout="BRP", C=C) + + # This is to know if plot a colorbar with optional labels + opt_cbar = nothing + if (showcpt) + if (cbar_label == "") opt_cbar = C + elseif (isa(cbar_label, String)) opt_cbar = (C, cbar_label) + elseif (isa(cbar_label, Tuple) && isa(cbar_label[1], String)) opt_cbar = (C, cbar_label...) + end + isa(cbar_label, String) && !(startswith(cbar_label, "xlabel=") || startswith(cbar_label, "ylabel=")) && + error("The labels in 'colorbar' option must start with 'xlabel=...' or 'ylabel=...'.") + isa(cbar_label, Tuple) && !((startswith(cbar_label[1], "xlabel=") || startswith(cbar_label[1], "ylabel=")) && + (startswith(cbar_label[2], "xlabel=") || startswith(cbar_label[2], "ylabel="))) && + @warn("colorbar labels MUST start with a `xlabel=` and/or `ylabel=` prefix.") + end + + if (!isempty(inset)) + Gs = interp_vslice(squeeze(slicecube(G, inset[2], x=[inset[1] G.range[2]], axis="y")), inc=interp) # inset South wall + Ge = interp_vslice(squeeze(slicecube(G, inset[1], y=[G.range[3] inset[2]], axis="x")), inc=interp) # inset East wall + pct_x = (G.range[2] - inset[1]) / (G.range[2] - G.range[1]) + pct_y = (inset[2] - G.range[3]) / (G.range[4] - G.range[3]) + resetGMT() # To remove "rememberings" left by the grdimage calls (they make annotations on West side??) + cubeplot(It, Is, Ie, R=@sprintf("%.12g/%.12g/%.12g/%.12g/%.12g/%.12g", G.range[1:4]..., G.range[7], G.range[8]), + hole=((grdimage(Gs, A=true, B=:none, layout="BRP", C=C), grdimage(Ge, A=true, B=:none, layout="BRP", C=C)), pct_x, pct_y), zsize=zsize, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, title=title, C=opt_cbar, show=show) + else + resetGMT() + cubeplot(It, Is, Ie, R=@sprintf("%.12g/%.12g/%.12g/%.12g/%.12g/%.12g", G.range[1:4]..., G.range[7], G.range[8]), + zsize=zsize, xlabel=xlabel, ylabel=ylabel, zlabel=zlabel, title=title, C=opt_cbar, show=show) + end +end + +# From https://stackoverflow.com/questions/54879412/get-the-mapping-from-each-element-of-input-to-the-bin-of-the-histogram-in-julia +# get the second output as in Matlab's histc +binindices(edges, data) = searchsortedlast.(Ref(edges), data) +# ------------------------------------------------------------------------------------------------- +""" + Gi = interp_vslice(G::GMTgrid; inc=0.0) -> GMTgrid + +Linearly interpolated the grid `G` along the "columns" (y coordinates). Normally `G` is a grid resulting +from a cube vertical slice where not uncomonly the resulting grid is not regular (cube layers are not equi-spaced). + +- `G`: A GMTgrid object with the grid to be interpolated +- `inc`: The interpolation increment. If == 0.0 and `G.y` is a constant spacing vector, nothing is done. + Otherwise, but still in the == 0.0 case, the grid interpolated with a pace equal to the minimum `G.y` spacing. + A `inc > 0.0` means the grid is interpolated at that increment. +""" +function interp_vslice(G::GMTgrid; inc=0.0) + # Do a linear interpolation along the "columns" (y coordinates) of the G grid + y = G.y + diff_y = diff(y) + if (inc == 0.0) + mi, ma = extrema(diff_y) + mi == ma && return G # No inc set and y has regular spacing, nothing to interpolate. + inc = min(mi, ma) + end + np = round(Int, (y[end] - y[1]) / inc + 1) + vk = linspace(y[1], y[end], np) # Vector of knots + ind = binindices(y, vk) # -> Vector{Int64} + append!(diff_y, 1) # Because we need an extra element at the end to make it same size as vk. That width of that extra bin is not used + f = (vk .- y[ind]) ./ diff_y[ind] + ff = (1.0 .- f) + mat = Matrix{eltype(G.z)}(undef, size(G,1), np) + @inbounds for k = 1:size(G,1) + @inbounds for n = 1:np-1 + mat[k,n] = G.z[k,ind[n]] * ff[n] + G.z[k,ind[n+1]] * f[n] + end + mat[k,np] = G.z[k,end] + end + mat2grid(mat, x=G.x, y=vk, layout=G.layout, is_transposed=true) +end From 624209960d5cddbc6d6f50725be21ff574833b96 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 10 Feb 2024 18:45:47 +0000 Subject: [PATCH 3/5] Add a -M alias --- src/pslegend.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pslegend.jl b/src/pslegend.jl index 554ed6e3b..a3130d3e8 100644 --- a/src/pslegend.jl +++ b/src/pslegend.jl @@ -20,7 +20,7 @@ Parameters - **F** | **box** :: [Type => Str | Number] `Arg=[+cclearances][+gfill][+i[[gap/]pen]][+p[pen]][+r[radius]][+s[[dx/dy/][shade]]]` Without further options, draws a rectangular border around the legend using *MAP_FRAME_PEN*. -- **M** :: [Type => Bool] +- **M** | **source** :: [Type => Bool] Modern mode only: - **S** | **scale** :: [Type => Number] @@ -55,7 +55,7 @@ function legend(cmd0::String="", arg1=nothing; first=true, kwargs...) opt_D = parse_type_anchor(d, "", [:D :pos :position], (map=("g", arg2str, 1), outside=("J", arg2str, 1), inside=("j", arg2str, 1), norm=("n", arg2str, 1), paper=("x", arg2str, 1), anchor=("", arg2str, 2), width=("+w", arg2str), justify="+j", spacing="+l", offset=("+o", arg2str)), 'j') - cmd = parse_these_opts(cmd, d, [[:C :clearance], [:M], [:S :scale], [:T :leg_file]]) + cmd = parse_these_opts(cmd, d, [[:C :clearance], [:M :source], [:S :scale], [:T :leg_file]]) SHOW_KWARGS[1] && legend_help() (opt_D == "") && error("The `position` argument is mandatory.") From f9327557bc7f9ede0f6637958857b36f8c402589 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 10 Feb 2024 18:47:11 +0000 Subject: [PATCH 4/5] Let title, xlabel, etc... be passed an empty string to mean no-op. --- src/common_options.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/common_options.jl b/src/common_options.jl index d97191949..fdc6b988a 100644 --- a/src/common_options.jl +++ b/src/common_options.jl @@ -991,11 +991,12 @@ function parse_B(d::Dict, cmd::String, opt_B__::String="", del::Bool=true)::Tupl # Let the :title and x|y_label be given on main kwarg list. Risky if used with NamedTuples way. t::String = "" # Use the trick to replace blanks by Char(127) (invisible) and undo it in extra_parse extra_par::String, tt::String, ep::String = "", "", "" - if (haskey(d, :title)) tt, extra_par = titlices(d, d[:title], title); t *= "+t" * tt end - if (haskey(d, :subtitle)) tt, ep = titlices(d, d[:subtitle], subtitle); t *= "+s" * tt; extra_par *= ep end - if (haskey(d, :xlabel)) tt, ep = titlices(d, d[:xlabel], xlabel); t *= " x+l" * tt; extra_par *= ep end - if (haskey(d, :ylabel)) tt, ep = titlices(d, d[:ylabel], ylabel); t *= " y+l" * tt; extra_par *= ep end - if (haskey(d, :zlabel)) tt, ep = titlices(d, d[:zlabel], zlabel); t *= " z+l" * tt; extra_par *= ep end + if (haskey(d, :title) && d[:title] != "") tt, extra_par = titlices(d, d[:title], title); t *= "+t" * tt end + if (haskey(d, :subtitle) && d[:subtitle] != "") tt, ep = titlices(d, d[:subtitle], subtitle);t *= "+s" * tt; extra_par *= ep end + if (haskey(d, :xlabel) && d[:xlabel] != "") tt, ep = titlices(d, d[:xlabel], xlabel); t *= " x+l" * tt; extra_par *= ep end + if (haskey(d, :ylabel) && d[:ylabel] != "") tt, ep = titlices(d, d[:ylabel], ylabel); t *= " y+l" * tt; extra_par *= ep end + if (haskey(d, :zlabel) && d[:zlabel] != "") tt, ep = titlices(d, d[:zlabel], zlabel); t *= " z+l" * tt; extra_par *= ep end + delete!(d, :title); delete!(d, :subtitle); delete!(d, :xlabel); delete!(d, :ylabel); delete!(d, :zlabel)# If == "" they were still there if (t != "") if (opt_B == "" && (val = find_in_dict(d, [:xaxis :yaxis :zaxis :xticks :yticks :zticks], false)[1] === nothing)) @@ -4703,7 +4704,7 @@ function interp_vec(x::AbstractVecOrMat{<:Real}, val::Real)::Float64 frac = (val - x[k]) / (x[k+1] - x[k]) return k + frac end -function interp_vec(x::AbstractVecOrMat{<:Real}, val::VecOrMat{<:Real}) +function interp_vec(x::AbstractVecOrMat{<:Real}, val::Union{VecOrMat{<:Real}, Tuple}) out = Vector{Float64}(undef, length(val)) for k = 1:numel(val) out[k] = interp_vec(x, val[k]) From 8f7a1b33edfcf6bfefc6cd1b25bae63b1d0f88b5 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 10 Feb 2024 18:47:44 +0000 Subject: [PATCH 5/5] Extend cube tests --- test/test_cube.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_cube.jl b/test/test_cube.jl index aedce4525..9b776d2b7 100644 --- a/test/test_cube.jl +++ b/test/test_cube.jl @@ -47,3 +47,8 @@ xyzw2cube("test_cube_ascii_colmaj.dat")[:,:,1] == [10.0 10.0 10.0; 10.0 10.0 xyzw2cube("test_cube_ascii_rowlevmaj.dat")[:,:,2] == [20.0 20.0 20.0; 20.0 20.0 20.0] xyzw2cube("test_cube_ascii_collevmaj.dat")[:,:,2] == [20.0 20.0 20.0; 20.0 20.0 20.0] xyzw2cube(gmtread("test_cube_ascii_rowmaj.dat"))[:,:,1] == [10.0 10.0 10.0; 10.0 10.0 10.0] + +cubeplot(cube, title="T annual", colorbar=("xlabel=bbb","ylabel=yy")) +cubeplot(cube, inset=(3,3)) +cubeplot(cube, top="@earth_relief_05m", inset=(3,3)) +cubeplot(cube, top="@earth_relief_05m", topshade=true)