From 43d983b785631c2d7a8961ab94984bda53761c5c Mon Sep 17 00:00:00 2001 From: Joaquim Date: Mon, 5 Feb 2024 22:36:59 +0000 Subject: [PATCH 1/3] Add sub2ind() Matlab clone function. --- src/utils.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index f046f4013..624d8f70a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -614,6 +614,11 @@ const numel = length dec2bin(n::Integer, mindigits::Int=0) = string(n, base=2, pad=mindigits) bin2dec(b::Union{AbstractString, Char}) = parse(Int, b, base=2) +function sub2ind(sz, args...) + linidx = LinearIndices(sz) + getindex.([linidx], args...) +end + function fileparts(fn::String) pato, ext = splitext(fn) pato, fname = splitdir(pato) From d50594525f418199acc5207dd2e07d9be760d3ca Mon Sep 17 00:00:00 2001 From: Joaquim Date: Mon, 5 Feb 2024 22:37:54 +0000 Subject: [PATCH 2/3] Improve & fixes on the cubeslice function. --- src/common_options.jl | 11 +++++- src/show_pretty_datasets.jl | 2 +- src/utils_types.jl | 79 +++++++++++++++++++++++-------------- test/test_misc.jl | 1 + test/test_new_projs.jl | 5 ++- 5 files changed, 65 insertions(+), 33 deletions(-) diff --git a/src/common_options.jl b/src/common_options.jl index ebe286174..d97191949 100644 --- a/src/common_options.jl +++ b/src/common_options.jl @@ -4692,17 +4692,24 @@ end # -------------------------------------------------------------------------------------------------- """ - frac = interp_vec(x, val) + frac = interp_vec(x, val) -> Float64 Returns the positional fraction that `val` ocupies in the `x` vector """ -function interp_vec(x, val) +function interp_vec(x::AbstractVecOrMat{<:Real}, val::Real)::Float64 (val < x[1] || val > x[end]) && error("Interpolating point ($val) is not inside the vector range [$(x[1]) $(x[end])].") k = 0 while(val < x[k+=1]) end frac = (val - x[k]) / (x[k+1] - x[k]) return k + frac end +function interp_vec(x::AbstractVecOrMat{<:Real}, val::VecOrMat{<:Real}) + out = Vector{Float64}(undef, length(val)) + for k = 1:numel(val) + out[k] = interp_vec(x, val[k]) + end + return out +end # -------------------------------------------------------------------------------------------------- function help_show_options(d::Dict) diff --git a/src/show_pretty_datasets.jl b/src/show_pretty_datasets.jl index 85625164e..18b7c70fd 100644 --- a/src/show_pretty_datasets.jl +++ b/src/show_pretty_datasets.jl @@ -253,7 +253,7 @@ For displaying data frame we do not want string representation of type to be longer than `maxwidth`. This function implements rules how type names are cropped if they are longer than `maxwidth`. """ -function compacttype(T::Type, maxwidth::Int) +function compacttype(@nospecialize(T::Type), maxwidth::Int) maxwidth = max(8, maxwidth) T === Any && return "Any" diff --git a/src/utils_types.jl b/src/utils_types.jl index 507bf3d98..efdb6c7fd 100644 --- a/src/utils_types.jl +++ b/src/utils_types.jl @@ -1238,7 +1238,7 @@ 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}}; axis="z") +function slicecube(G::GMTgrid, slice::Union{Int, AbstractVector{<:Int}}; x::VecOrMat{<:Real}=Float64[], y::VecOrMat{<:Real}=Float64[], axis="z") # Method that slices grid cubes. SLICE is the row|col|layer number. AXIS picks the axis to be sliced (ndims(G) < 3 || size(G,3) < 2) && error("This is not a cube grid.") _axis = lowercase(string(axis)) @@ -1251,57 +1251,66 @@ function slicecube(G::GMTgrid, slice::Union{Int, AbstractVector{<:Int}}; axis="z (!isvec && slice[end] > this_size) && error("Last slice number ($(slice[end])) is larger than grid size ($this_size)") 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=(G.layout[2] == 'R')) + G_ = mat2grid(G[:,:,slice], G.x, G.y, 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") - if (G.layout[2] == 'C') G_ = mat2grid(G[slice,:,:], G.v, G.x, reg=G.registration, names=G.names) - else G_ = mat2grid(G[:,slice,:], G.x, G.v, reg=G.registration, is_transposed=true, names=G.names) + 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) end - G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = "", "", 0, 0, G.layout - G_.v = G_.y; G_.y = isvec ? G.y[slice] : [G.y[slice]]; # Shift coords vectors since mat2grid doesn't know how-to. + G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = "", "", 0, 0, "TRB" else - if (G.layout[2] == 'C') G_ = mat2grid(G[:,slice,:], G.v, G.y, reg=G.registration, names=G.names) - else G_ = mat2grid(G[slice,:,:], G.y, G.v, reg=G.registration, is_transposed=true, names=G.names) + 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) end - G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = "", "", 0, 0, G.layout - G_.v = G_.y; G_.y = G_.x; G_.x = isvec ? G.x[slice] : [G.x[slice]]; + G_.proj4, G_.wkt, G_.epsg, G_.geog, G_.layout = "", "", 0, 0, "TRB" end return G_ end -function slicecube(G::GMTgrid, slice::AbstractFloat; axis="z") +function slicecube(G::GMTgrid, slice::AbstractFloat; x::VecOrMat{<:Real}=Float64[], y::VecOrMat{<:Real}=Float64[], axis="z") # Method that slices grid cubes. SLICE is the x|y|z coordinate where to slice. AXIS picks the axis to be sliced (ndims(G) < 3 || size(G,3) < 2) && error("This is not a cube grid.") _axis = lowercase(string(axis)) which_coord_vec = (_axis == "z") ? G.v : (_axis == "y" ? G.y : G.x) - isempty(which_coord_vec) && (which_coord_vec = 1.0:size(G,3)) # To at least have something. - x = interp_vec(which_coord_vec, slice) - layer = trunc(Int, x) - frac = x - layer # Used for layer interpolation. - (frac < 0.1) && return slicecube(G, layer) # If 'slice' is within 10% of lower or upper layer just - (frac > 0.9) && return slicecube(G, layer+1) # return that layer and do no interpolation between layers. + isempty(which_coord_vec) && (which_coord_vec = collect(1.0:size(G,3))) # To at least have something. + xf = interp_vec(which_coord_vec, slice) + layer = trunc.(Int, xf) + frac = xf .- layer # Used for layer interpolation. + 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) + l = layer + colmajor = (G.layout[2] == 'C') if (_axis == "z") 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=(G.layout[2] == 'R')) + 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") - if (G.layout[2] == 'C') mat = G[layer,:,:] .+ (G[layer+1,:,:] .- G[layer,:,:]) .* frac - else mat = G[:,layer,:] .+ (G[:,layer+1,:] .- G[:,layer,:]) .* frac # from GDAL + 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) 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, G.v, reg=G.registration, is_transposed=(G.layout[2] == 'R')) - G_.v = G_.y; G_.y = [Float64(slice)] # Shift coords vectors since mat2grid doesn't know how-to. + G_ = mat2grid(mat, G.x[ix], [G.y[l]], G.v, reg=G.registration, is_transposed=!colmajor, layout="TRB") else - if (G.layout[2] == 'C') mat = G[:,layer,:] .+ (G[:,layer+1,:] .- G[:,layer,:]) .* frac - else mat = G[layer,:,:] .+ (G[layer+1,:,:] .- G[layer,:,:]) .* frac # from GDAL + 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) 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.y, G.v, reg=G.registration, is_transposed=(G.layout[2] == 'R')) - G_.v = G_.y; G_.y = G_.x; G_.x = [Float64(slice)] + G_ = mat2grid(mat, [G.x[l]], G.y[iy], G.v, reg=G.registration, is_transposed=!colmajor, layout="TRB") end - G_.layout = G.layout return G_ end @@ -1353,6 +1362,18 @@ function slicecube(GI::GItype; slice::Int=0, α=0.0, angle=0.0, axis="x", cmap=G end const cubeslice = slicecube # I'm incapable of remembering which one it is. + +# --------------------------------------------------------------------------------------------------- +""" +""" +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) + end +end + # --------------------------------------------------------------------------------------------------- """ xyzw2cube(fname::AbstractString; zcol::Int=4, datatype::DataType=Float32, proj4::String="", wkt::String="", @@ -1865,8 +1886,8 @@ function grdimg_hdr_xy(mat, reg, hdr, x=Float64[], y=Float64[], is_transposed=fa x = collect(range(x[1], stop=x[2], length=nx+reg)) y = collect(range(y[1], stop=y[2], length=ny+reg)) end - x_inc = (x[end] - x[1]) / (nx - one_or_zero) - y_inc = (y[end] - y[1]) / (ny - one_or_zero) + x_inc = (x[end] - x[1]) / (nx - one_or_zero); isnan(x_inc) && (x_inc = 0.0) # When vertical slices + y_inc = (y[end] - y[1]) / (ny - one_or_zero); isnan(y_inc) && (y_inc = 0.0) zmin::Float64, zmax::Float64 = extrema_nan(mat) hdr = Float64.([x[1], x[end], y[1], y[end], zmin, zmax]) elseif (isempty(hdr)) diff --git a/test/test_misc.jl b/test/test_misc.jl index 1e0646bd0..4ca0c2aa0 100644 --- a/test/test_misc.jl +++ b/test/test_misc.jl @@ -51,6 +51,7 @@ GMT.resetGMT() @test GMT.bin2dec("10111") == 23 @test GMT.dec2bin(23) == "10111" + @test GMT.sub2ind((3,3), [1 2 3 1], [2 2 2 3]) == [4 5 6 7] A = [3 1; 3 3; 1 3; 3 2; 2 3; 1 1; 1 2; 2 3; 3 3; 3 3]; C, ia, ic = GMT.uniq(A; dims=1); diff --git a/test/test_new_projs.jl b/test/test_new_projs.jl index 95708a712..54711030d 100644 --- a/test/test_new_projs.jl +++ b/test/test_new_projs.jl @@ -16,7 +16,10 @@ end cubeplot(mat2img(rand(UInt8,8,8,3)), p=(145,30), front=false, zsize=10, notop=true, back=true) - cubeplot(mat2img(rand(UInt8,8,8,3)), mat2img(rand(UInt8,8,8,3))) + I1 = mat2img(rand(UInt8,8,8,3)); + I2 = mat2img(rand(UInt8,8,8,3)); + I3 = mat2img(rand(UInt8,8,8,3)); + cubeplot(I1, I2, I3, inset=((I1,I2), .5) cl = coastlinesproj(proj="fouc"); grid = graticules(cl); From 6274284d239b9f1311a5f392734dd3833b955a68 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Mon, 5 Feb 2024 22:50:35 +0000 Subject: [PATCH 3/3] Fix syntax error --- test/test_new_projs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_new_projs.jl b/test/test_new_projs.jl index 54711030d..381f78940 100644 --- a/test/test_new_projs.jl +++ b/test/test_new_projs.jl @@ -19,7 +19,7 @@ I1 = mat2img(rand(UInt8,8,8,3)); I2 = mat2img(rand(UInt8,8,8,3)); I3 = mat2img(rand(UInt8,8,8,3)); - cubeplot(I1, I2, I3, inset=((I1,I2), .5) + cubeplot(I1, I2, I3, inset=((I1,I2), .5)) cl = coastlinesproj(proj="fouc"); grid = graticules(cl);