diff --git a/src/GMT.jl b/src/GMT.jl index cbb646709..75eb9b2cf 100644 --- a/src/GMT.jl +++ b/src/GMT.jl @@ -204,6 +204,7 @@ include("utils_types.jl") include("grd_operations.jl") include("common_options.jl") const global LEGEND_TYPE = [legend_bag()] # To store Legends info +include("circfit.jl") include("custom_symb_funs.jl") include("gmtbegin.jl") include("blendimg.jl") @@ -387,7 +388,7 @@ using GMT.Laszip resetGMT() end -Base.precompile(Tuple{typeof(upGMT),Bool, Bool}) # Here it doesn't print anything. +#Base.precompile(Tuple{typeof(upGMT),Bool, Bool}) # Here it doesn't print anything. function __init__(test::Bool=false) if !isfile(libgmt) || (!Sys.iswindows() && (!isfile(libgdal) || !isfile(libproj))) diff --git a/src/blendimg.jl b/src/blendimg.jl index ea472230c..5d877db12 100644 --- a/src/blendimg.jl +++ b/src/blendimg.jl @@ -46,7 +46,7 @@ function blendimg!(color::GMTimage, shade::GMTimage; mode="", transparency=0.5, end # --------------------------------------------------------------------------------------------------- -function blend_gcalc(color::GMTimage{UInt8, 3}, shade::GMTimage{UInt8, 2}; new=false) +function blend_gcalc(color::GMTimage{UInt8, 3}, shade::GMTimage{UInt8, 2}; new=false)::GMTimage blend = (new) ? Array{UInt8,3}(undef, size(shade,1), size(shade,2), 3) : color.image @@ -79,7 +79,7 @@ function blend_gcalc(color::GMTimage{UInt8, 3}, shade::GMTimage{UInt8, 2}; new=f end # --------------------------------------------------------------------------------------------------- -function blend_blend(img1::GMTimage, img2::GMTimage; new=false, transparency=0.5) +function blend_blend(img1::GMTimage, img2::GMTimage; new=false, transparency=0.5)::GMTimage # This method blends two UInt8 images with transparency @assert eltype(img1) == eltype(img2) @assert length(img1) == length(img2) @@ -105,7 +105,7 @@ function blend_blend(img1::GMTimage, img2::GMTimage; new=false, transparency=0.5 end # --------------------------------------------------------------------------------------------------- -function blend_PS(A::GMTimage{UInt8, 3}, B::GMTimage; mode="LinearBurn", new=false, white=220, burn_up=180, screen_low=220) +function blend_PS(A::GMTimage{UInt8, 3}, B::GMTimage; mode="LinearBurn", new=false, white=220, burn_up=180, screen_low=220)::GMTimage # PhotoShop blend modes (a few) blend = (new) ? Array{UInt8,size(A,3)}(undef, size(A)) : A.image @@ -201,7 +201,7 @@ Optionally set also `contrast` and/or `brightness` ### Returns A GMT intensity Image """ -function gammacorrection(I::GMTimage, gamma; contrast=[0.0, 1.0], brightness=[0.0, 1.0]) +function gammacorrection(I::GMTimage, gamma; contrast=[0.0, 1.0], brightness=[0.0, 1.0])::GMTimage @assert 0.0 <= contrast[1] < 1.0; @assert 0.0 < contrast[2] <= 1.0; @assert contrast[2] > contrast[1] @assert 0.0 <= brightness[1] < 1.0; @assert 0.0 < brightness[2] <= 1.0; @assert brightness[2] > brightness[1] @@ -255,7 +255,7 @@ http://www.textureshading.com/Home.html ### Returns A UInt8 (or 16) GMT Image """ -function texture_img(G::GMTgrid; detail=1.0, contrast=2.0, uint16=false, intensity=false) +function texture_img(G::GMTgrid; detail=1.0, contrast=2.0, uint16=false, intensity=false)::GMTimage # Here we have a similar problem with the memory layout as described in gmt2ds(). Specialy with the # variants of the TRB layout. Other than that, it's strange (probably to dive in the C code to relearn) # why the memory layout of the BCB mode needs to changed in the way we do below. diff --git a/src/common_options.jl b/src/common_options.jl index d278b2009..8b120553f 100644 --- a/src/common_options.jl +++ b/src/common_options.jl @@ -3954,6 +3954,18 @@ function isFV(D)::Bool (D[1].geom == wkbPoint || D[1].geom == wkbPointZ) && eltype(D[2].data) <: Integer end +# --------------------------------------------------------------------------------------------------- +""" + is_gridtri(D)::Bool + +Check if D is a Vector{GMTdataset} produced by the gridtri() function. +""" +function is_gridtri(D)::Bool + (!isa(D, Vector{<:GMTdataset}) || isempty(D[1].comment)) && return false + (!contains(D[1].comment[1], "vwall") && !contains(D[1].comment[1], "gridtri")) && return false + return true +end + # --------------------------------------------------------------------------------------------------- function find_data(d::Dict, cmd0::String, cmd::String, args...) # ... diff --git a/src/gmt_main.jl b/src/gmt_main.jl index 930d36694..33ab67728 100644 --- a/src/gmt_main.jl +++ b/src/gmt_main.jl @@ -1180,8 +1180,8 @@ function palette_init(API::Ptr{Nothing}, cpt::GMTcpt)::Ptr{GMT_PALETTE} n_colors = n_colors - 1; # Number of CPT slices end - P = convert(Ptr{GMT_PALETTE}, GMT_Create_Data(API, GMT_IS_PALETTE, GMT_IS_NONE, 0, pointer([n_colors]), - NULL, NULL, 0, 0, NULL)) + ptr = Ref{UInt64}(n_colors) + P = convert(Ptr{GMT_PALETTE}, GMT_Create_Data(API, GMT_IS_PALETTE, GMT_IS_NONE, 0, ptr, NULL, NULL, 0, 0, NULL)) (n_colors > 100000) && @warn("Que exagero de cores") # Just to protect n_colors to be GC'ed before here Pb::GMT_PALETTE = unsafe_load(P) # We now have a GMT_PALETTE @@ -1251,7 +1251,7 @@ function ps_init(API::Ptr{Nothing}, ps, dir::Integer)::Ptr{GMT_POSTSCRIPT} (!isa(ps, GMTps)) && error("Expected a PS structure for input") # Passing dim[0] = 0 since we dont want any allocation of a PS string - pdim = pointer([0]) + pdim = Ref{UInt64}(0) P = convert(Ptr{GMT_POSTSCRIPT}, GMT_Create_Data(API, GMT_IS_POSTSCRIPT, GMT_IS_NONE, 0, pdim, NULL, NULL, 0, 0, NULL)) P0::GMT_POSTSCRIPT = unsafe_load(P) # GMT_POSTSCRIPT diff --git a/src/imshow.jl b/src/imshow.jl index 38c05744b..1feebece2 100644 --- a/src/imshow.jl +++ b/src/imshow.jl @@ -40,11 +40,11 @@ function imshow(arg1, x::AbstractVector{Float64}=Float64[], y::AbstractVector{Fl see = ((val = find_in_kwargs(kw, [:show])[2]) === nothing) ? true : (val != 0) # No explicit 'show' keyword means show=true function isplot3(kw) - call_plot3 = false + _call_plot3 = false opt_p = find_in_kwargs(kw, [:p :view :perspective])[1] - (isa(opt_p, String) && contains(opt_p, '/')) && (call_plot3 = true) - ((isa(opt_p, Tuple) || isa(opt_p, VMr)) && length(opt_p) > 1) && (call_plot3 = true) - return call_plot3 + (isa(opt_p, String) && contains(opt_p, '/')) && (_call_plot3 = true) + ((isa(opt_p, Tuple) || isa(opt_p, VMr)) && length(opt_p) > 1) && (_call_plot3 = true) + return _call_plot3 end is_image, call_img, call_grd = false, false, false diff --git a/src/mapproject.jl b/src/mapproject.jl index 079be9264..aa00f560a 100644 --- a/src/mapproject.jl +++ b/src/mapproject.jl @@ -101,12 +101,16 @@ function mapproject(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...) (dbg_print_cmd(d, cmd) !== nothing) && return cmd gmt("mapproject " * cmd) else - #cmd, got_fname, arg1, arg2 = find_data(d, cmd0, cmd, arg1, arg2) - (isa(arg1, Vector{<:Real})) && (arg1 = reshape(arg1, 1, :)) # User may have send a vector but GMT expects matrix + #cmd, got_fname, arg1, arg2 = find_data(d, cmd0, cmd, arg1, arg2) + (isa(arg1, Vector{<:Real})) && (arg1 = reshape(arg1, 1, :)) # User may have send a vector but GMT expects matrix R = common_grd(d, cmd0, cmd, "mapproject ", arg1, arg2) contains(cmd, " -I") && (isa(R, GMTdataset) ? (R.proj4 = prj4WGS84) : (R[1].proj4 = prj4WGS84)) - contains(opt_J, "+proj") && (isa(R, GMTdataset) ? (R.proj4 = opt_J[4:end]) : (R[1].proj4 = opt_J[4:end])) - R + if contains(opt_J, "+proj") + prj = opt_J[4:end] + contains(prj, '+') && !contains(prj, ' ') && (prj = replace(prj, "+" => " +")) # Can't store a 'glued' -J + (isa(R, GMTdataset) ? (R.proj4 = prj) : (R[1].proj4 = prj)) + end + R end end diff --git a/src/psxy.jl b/src/psxy.jl index 921694f3b..e784c7ba6 100644 --- a/src/psxy.jl +++ b/src/psxy.jl @@ -58,6 +58,9 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D:: if (is3D && isFV(arg1)) # case of 3D faces arg1 = deal_faceverts(arg1, d, opt_p) !haskey(d, :aspect3) && (d[:aspect3] = "equal") # Needs thinking + elseif is_gridtri + arg1 = sort_visible_triangles(arg1) + is_in_dict(d, [:Z :level :levels]) === nothing && (d[:Z] = tri_z(arg1)) end isa(arg1, GMTdataset) && (arg1 = with_xyvar(d, arg1)) # See if we have a column request based on column names @@ -371,7 +374,7 @@ function deal_faceverts(arg1, d, opt_p) elev = length(spl) > 1 ? parse(Float64, spl[2]) : 90.0 arg1, dotprod = visible_faces(arg1, [sind(az) * cosd(elev), cosd(az) * cosd(elev), sind(elev)]) if (is_in_dict(d, [:G :fill]) === nothing) # If fill not set we use the dotprod and a gray CPT to set the fill - d[:Z] = dotprod + is_in_dict(d, [:Z :level :levels]) === nothing && (d[:Z] = dotprod) (is_in_dict(d, CPTaliases) === nothing) && (d[:C] = gmt("makecpt -T0/1 -C150,210")) # Users may still set a CPT end return arg1 @@ -380,11 +383,12 @@ end # --------------------------------------------------------------------------------------------------- function deal_gridtri!(arg1, d)::Bool # Deal with the situation where we are plotting triangulated grids made by grid2tri() - ((!isa(arg1, Vector{<:GMTdataset}) || isempty(arg1[1].comment) || (arg1[1].comment[1] != "vwall" && arg1[1].comment[1] != "gridtri"))) && - return false + !is_gridtri(arg1) && return false is_in_dict(d, [:G :fill]) === nothing && (d[:G] = "+z") if (is_in_dict(d, CPTaliases) === nothing) - C = gmt("makecpt -T$(arg1[1].ds_bbox[5]-1e-8)/$(arg1[1].ds_bbox[6]+1e-8)/+n255 -Cturbo") + opt_T = (contains(arg1[1].comment[1], "vwall+gridtri_top")) ? "-T$(arg1[1].comment[2])/$(arg1[1].ds_bbox[6]+1e-8)/+n255" : + "-T$(arg1[1].ds_bbox[5]-1e-8)/$(arg1[1].ds_bbox[6]+1e-8)/+n255" + C = gmt("makecpt -Cturbo " * opt_T) C.bfn[2, :] .= 0.7 # Set the foreground color used by the vertical wall d[:C] = C end @@ -1426,3 +1430,60 @@ function visible_faces(FV::Vector{<:GMTdataset}, view_vec) F2 = mat2ds(F.data[is_vis, :]) return [V,F2], proj[is_vis] end + +# --------------------------------------------------------------------------------------------------- +function sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfact=1.0)::Vector{GMTdataset} + azim, elev = parse.(Float64, split(CURRENT_VIEW[1][4:end], '/')) + sin_az, cos_az, sin_el = sind(azim), cosd(azim), sind(elev) + prj, wkt, epsg = Dv[1].proj4, Dv[1].wkt, Dv[1].epsg + top_comment = Dv[1].comment # save this for later restore as it can be used by tri_z() to detect vertical walls + + (del_hidden != 1 && contains(Dv[1].comment[1], "vwall")) && (del_hidden = true) # If have vwalls, need to del invis + if (del_hidden == 1) # Remove the triangles that are not visible from the normal view_vec + t = isgeog(Dv) ? mapproject(Dv, J="t$((Dv[1].ds_bbox[1] + Dv[1].ds_bbox[2])/2)/1:1", C=true, F=true) : Dv + view_vec = [sin_az * cosd(elev), cos_az * cosd(elev), sin_el] + #zs_t = CTRL.pocket_J[3][5:end] # Ex: CTRL.pocket_J[3] = " -JZ3" + is_vis = [dot(facenorm(t[k].data, zfact=zfact, normalize=false), view_vec) > 0 for k in eachindex(t)] + #@show(length(t),sum(is_vis)) + Dv = Dv[is_vis] + end + + # ---------------------- Now sort by distance to the viewer ---------------------- + Dc = gmtspatial(Dv, Q=true, o="0,1") + dists = [(Dc.data[1,1] * sin_az + Dc.data[1,2] * cos_az, (Dv[1].bbox[5] + Dv[1].bbox[6]) / 2 * sin_el)] + for k = 2:size(Dc, 1) + push!(dists, (Dc.data[k,1] * sin_az + Dc.data[k,2] * cos_az, (Dv[k].bbox[5] + Dv[k].bbox[6]) / 2 * sin_el)) + end + + ind = sortperm(dists) # Sort in growing distances. + Dv = Dv[ind] + set_dsBB!(Dv) + Dv[1].proj4 = prj; Dv[1].wkt = wkt; Dv[1].epsg = epsg # Because first triangle may have been deleted or reordered. + Dv[1].comment = top_comment # Restore the original comment that holds info abot this gridtri dataset + return Dv +end + +# --------------------------------------------------------------------------------------------------- +function tri_normals(Dv::Vector{<:GMTdataset}; zfact=1.0) + t = isgeog(Dv) ? mapproject(Dv, J="t$((Dv[1].ds_bbox[1] + Dv[1].ds_bbox[2])/2)/1:1", C=true, F=true) : Dv + Dc = gmtspatial(Dv, Q=true) + nx = Vector{Float64}(undef, length(t)); ny = copy(nx); nz = copy(nx) + for k in eachindex(t) + nx[k], ny[k], nz[k] = facenorm(t[k].data; zfact=zfact) + Dc.data[k,3] = (t[k].bbox[5] + t[k].bbox[6]) / 2 # Reuse the area column to store the triangle mean height + end + return Dc, nx, ny, nz +end + +# --------------------------------------------------------------------------------------------------- +function quiver3(Dv::Vector{<:GMTdataset}; first=true, zfact=1.0, kwargs...) + Dc, nx, ny, nz = tri_normals(Dv, zfact=zfact) + mat = fill(NaN, size(Dc,1) * 3 - 1, 3) + for k in eachindex(Dv) + kk = (k-1) * 3 + mat[kk+=1,1], mat[kk,2], mat[kk,3] = Dc.data[k,1], Dc.data[k,2], Dc.data[k,3] + mat[kk+=1,1], mat[kk,2], mat[kk,3] = Dc.data[k,1]+nx[k]/2, Dc.data[k,2]+ny[k]/2, Dc.data[k,3]+nz[k]/2 + end + D = mat2ds(mat, geom=wkbLineStringZ) + common_plot_xyz("", D, "plot3d", first, true, kwargs...) +end diff --git a/src/triangulate.jl b/src/triangulate.jl index 3e434562e..d59300dcc 100644 --- a/src/triangulate.jl +++ b/src/triangulate.jl @@ -155,30 +155,37 @@ Since this is a `plot3d` _avatar_ all options in this function are those of the x,y,z = GMT.peaks(N=45, grid=false); trisurf([x[:] y[:] z[:]], pen=0.5, show=true) """ -function trisurf(in::Union{Matrix, GDtype}; gdal=true, first::Bool=true, kw...) - (size(in, 2) < 3) && error("'trisurf' needs at least 3 columns in input") +function trisurf(in::Union{Matrix, GDtype}; first::Bool=true, gdal::Bool=false, kw...) + # Keep the 'gdal' kwarg for backward compatibility (no longer used) d = KW(kw) first && (d[:aspect] = get(d, :aspect, "equal")) - d[:p] = get(d, :p, "135/30") - D = gdal ? delaunay(in, 0.0, false) : triangulate(in, S=true, Z=true) - D, Zs = sort_by_depth(D) # Sort the triangles by depth. Deepest first. Otherwise the plot is always a mess - d[:Z] = Zs + (!first && is_in_dict(d, [:p :projection]) === nothing ) && (d[:p] = "217.5/30")# This case is not tested in psxy + is_gridtri(in) && return common_plot_xyz("", in, "plot3d", first, true, d...) # Also works for gridtri's + + (size(in, 2) < 3) && error("This 'trisurf' call needs at least 3 columns in input") + D = triangulate(in, S=true, Z=true) + D[1].comment = ["gridtri_top"] common_plot_xyz("", D, "plot3d", first, true, d...) end -trisurf!(in::Union{Matrix, GDtype}; gdal=true, kw...) = trisurf(in; gdal=gdal, first=false, kw...) +trisurf!(in::Union{Matrix, GDtype}; kw...) = trisurf(in; first=false, kw...) -# --------------------------------------------------------------------------------------------------- -function sort_by_depth(D::Vector{<:GMTdataset}) - # Sort the triangulation `D` by depth (Z). Deepest first. - Zs = Vector{Float64}(undef, length(D)) - for k = 1:numel(Zs) - Zs[k] = (D[k].bbox[5] + D[k].bbox[6]) / 2 - end - ind = sortperm(Zs) # Sort in groing z. Needed to sort the triangles too otherwise is a mess. - ds_bbox = D[1].ds_bbox - D = D[ind]; D[1].ds_bbox = ds_bbox - return D, Zs[ind] +""" + trisurf(G, G2=nothing; bottom=false, downsample=0, level=false, ratio=0.01, + thickness=0.0, wall_only=false, top_only=false, geog=false, kw...) + +Short form for the sequence ``D = grid2tri(...)`` followed by ``viz(D)``. See the documentation of +``grid2tri`` for more details. +""" +function trisurf(G::Union{GMTgrid, String}, G2=nothing; first::Bool=true, thickness=0.0, level=false, downsample=0, + ratio=0.01, bottom=false, wall_only=false, top_only=false, geog=false, kw...) + D = grid2tri(G, G2, thickness=thickness, level=level, downsample=downsample, ratio=ratio, bottom=bottom, + wall_only=wall_only, top_only=top_only, geog=geog) + trisurf(D; first=first, kw...) end +trisurf!(G::Union{GMTgrid, String}, G2=nothing; thickness=0.0, level=false, downsample=0, ratio=0.01, + bottom=false, wall_only=false, top_only=false, geog=false, kw...) = + trisurf(G, G2; first=false, thickness=thickness, level=level, downsample=downsample, ratio=ratio, + bottom=bottom, wall_only=wall_only, top_only=top_only, geog=geog, kw...) # --------------------------------------------------------------------------------------------------- # D = grid2tri("cam_slab2_dep_02.24.18.grd", "cam_slab2_thk_02.24.18.grd") @@ -199,13 +206,13 @@ NOTE: The `G` grid should have a _out-skirt_ of NaNs, otherwise just use ``grdvi ### Parameters - `G`: A GMTgrid object or a grid file name representing the surface to be triangulated. -- `G2`: An optional second grid (ot file name) representing the bottom surface of the layer. Using this +- `G2`: An optional second grid (or file name) representing the bottom surface of the layer. Using this option makes the `thickness` option be ignored. ### Keywords - `bottom`: If true, fully close the body with the bottom surface. By default we don't do this because that - surface is not visible whem elevation view angle is positive. But we may want this if later we want to save - this mesh in STL for importing in a 3D viewer software. + surface is often not visible whem elevation view angle is positive. But we may want this if later we want + to save this mesh in STL for importing in a 3D viewer software. - `layer`: If true, we interpret `thickness` option as meaning a contant level value. That is, the vertical wall is computed from the sides of `G` and a constant level provided via the `thickness` option. @@ -234,47 +241,44 @@ NOTE: The `G` grid should have a _out-skirt_ of NaNs, otherwise just use ``grdvi ### Returns A vector of GMTdataset with the triangulated surface and vertical wall, or just the wall or the full closed body. """ -function grid2tri(G, G2=nothing; thickness=0.0, level=false, downsample=0, ratio=0.01, bottom=false, +function grid2tri(G::Union{GMTgrid, String}, G2=nothing; thickness=0.0, level=false, downsample=0, ratio=0.01, bottom=false, wall_only=false, top_only=false, geog=false) - (!isa(G2, GMTgrid) && !isa(G2, String) && thickness <= 0.0 && top_only == 0) && - error("Must provide a bottom grid or a thickness for this layer.") + (!isa(G2, GMTgrid) && !isa(G2, String) && thickness <= 0.0 && top_only == 0) && (top_only = true) (wall_only != 0) && (bottom = false) (top_only != 0) && (G2 = nothing; bottom = false; wall_only = false) Dbnd_t, Dpts = gridhull(G; downsample=downsample, ratio=ratio) # Compute the top concave hull - if (wall_only == 0) # (wall_only = 0 means we have to compute the top surface) - Dt_t = triangulate(Dpts, S="+za", Z=true) # Triangulation of top surface + if (wall_only == 0) # (wall_only == 0 means we have to compute the top surface) + Dt_t = triangulate(Dpts, S=true, Z=true) # Triangulation of top surface Dc = gmtspatial(Dt_t, Q=true, o="0,1") # Compute the polygon centroids ind = (Dc in Dbnd_t) .== 1 Dt_t = Dt_t[ind] # Delete the triangles outside the concave hull - Dt_t = sort_by_depth(Dt_t)[1] # Sort the triangles by depth. Deepest first. Dt_t[1].ds_bbox = Dt_t[1].bbox # Because we may have deleted first Dt_t Dt_t[1].proj4 = Dbnd_t.proj4 Dt_t[1].geom = wkbPolygonZM - Dt_t[1].comment = ["gridtri"] # To help recognize this type of DS + Dt_t[1].comment = ["gridtri_top"] # To help recognize this type of DS set_dsBB!(Dt_t, false) end - if (isa(G2, GMTgrid) || isa(G2, String)) + if (isa(G2, GMTgrid) || isa(G2, String)) # Gave a bottom surface Dbnd_b, Dpts = gridhull(G2; downsample=downsample, ratio=ratio) # Compute the bottom concave hull if (bottom == 1) - Dt_b = triangulate(Dpts, S="+za", Z=true) # Triangulation of bottom surface + Dt_b = triangulate(Dpts, S=true, Z=true) # Triangulation of bottom surface Dt_b = Dt_b[ind] # Delete the triangles outside the concave hull (reuse the same 'ind') - Dt_b = sort_by_depth(Dt_b)[1] # Sort the triangles by depth. Deepest first. Dt_b[1].ds_bbox = Dt_b[1].bbox # Because we may have deleted first Dt Dt_b[1].geom = wkbPolygonZM - Dt_b[1].comment = ["gridtri"] # To help recognize this type of DS + Dt_b[1].comment = ["gridtri_bot"] # To help recognize this type of DS for k = 1:numel(Dt_b) # Must reverse the order of the vertices so normals point outward Dt_b[k][3,1], Dt_b[k][2,1] = Dt_b[k][2,1], Dt_b[k][3,1] Dt_b[k][3,2], Dt_b[k][2,2] = Dt_b[k][2,2], Dt_b[k][3,2] Dt_b[k][3,3], Dt_b[k][2,3] = Dt_b[k][2,3], Dt_b[k][3,3] end - set_dsBB!(Dt_b, false) end Dwall = vwall(Dbnd_t, view(Dbnd_b, :, 3)) if (bottom == 1) append!(Dt_b, Dwall, Dt_t) # If including bottom too, start with it and add the wall and top. - (geog == 0) ? refsystem_A2B!(Dbnd_t, Dt_b) : (Dt_b[1].proj4 = prj4WGS84) # Set ref sys + (geog == 0) ? refsystem_A2B!(Dbnd_t, Dt_b) : (Dt_b[1].proj4 = prj4WGS84) # Set ref sys + set_dsBB!(Dt_b, false) return Dt_b end else @@ -285,8 +289,13 @@ function grid2tri(G, G2=nothing; thickness=0.0, level=false, downsample=0, ratio Dwall = vwall(Dbnd_t, thickness, level != 0) end - (wall_only == 0) && append!(Dwall, Dt_t) + if (wall_only == 0) # That is vertical wall + top surface + append!(Dwall, Dt_t) + Dwall[1].comment[1] = "vwall+gridtri_top" + append!(Dwall[1].comment, ["$(Dt_t[1].ds_bbox[5])"]) # Save the top surface minimum to help making a default cpt. + end (geog == 0) ? refsystem_A2B!(Dbnd_t, Dwall) : (Dwall[1].proj4 = prj4WGS84) # Set ref sys + set_dsBB!(Dwall, false) return Dwall end @@ -300,8 +309,8 @@ end - `downsample`: Downsample the input grid by `downsample` times. - `ratio`: The ratio of the concave hull to the convex hull. """ -function gridhull(G; downsample::Int=0, ratio=0.01) - _G = isa(G, String) ? gmtread(G) : G +function gridhull(G::Union{GMTgrid, String}; downsample::Int=0, ratio=0.01) + _G::GMTgrid = isa(G, String) ? gmtread(G) : G prj4, wkt, epsg = _G.proj4, _G.wkt, _G.epsg # Save the original projection because currently grdsample looses them (downsample > 1) && (_G = grdsample(_G, I="$(div(size(_G.z,2),downsample)+1)+n/$(div(size(_G.z,1),downsample)+1)+n", V="q")) V = grd2xyz(_G, s=true) # Convert to x,y,z while dropping the NaNs @@ -350,7 +359,7 @@ function vwall(Bt::Union{Matrix{<:Real}, GMTdataset}, Bb::Union{Matrix{<:Real}, Bb[k,1] Bb[k,2] Bb[k,3]; Bt[k,1] Bt[k,2] Bt[k,3]]) Twall[kk+1] = GMTdataset(data=[Bt[k+1,1] Bt[k+1,2] Bt[k+1,3]; Bb[k+1,1] Bb[k+1,2] Bb[k+1,3]; Bb[k,1] Bb[k,2] Bb[k,3]; Bt[k+1,1] Bt[k+1,2] Bt[k+1,3]]) - Twall[kk].header = Twall[kk+1].header = "-Z$(Bt[k,3]+1e6)" # Can be used, e.g., with the foreground color of a CPT + Twall[kk].header = Twall[kk+1].header = "W " # To allow identifying this as a vertical wall end set_dsBB!(Twall) isa(Bt, GMTdataset) && (Twall[1].proj4 = Bt.proj4) @@ -358,3 +367,28 @@ function vwall(Bt::Union{Matrix{<:Real}, GMTdataset}, Bb::Union{Matrix{<:Real}, Twall[1].comment = ["vwall"] # To help recognize this as a vertical wall. return Twall end + +# --------------------------------------------------------------------------------------------------- +""" + Z = tri_z(D::Vector{<:GMTdataset}) + +Get the half of the elevation for each 3D polygon in the vector of `D`. Note: this is NOT the average +of vertices elevation, but the elevation at the midpoint of the polygon. + +In case the `D` is vector contains also a vertical wall, signaled by the comment starting with "vwall", +we set the elevation to 1e6 so that we can send these triangles to the foreground color of a CPT. +""" +function tri_z(D::Vector{<:GMTdataset})::Vector{Float64} + Zs = Vector{Float64}(undef, length(D)) + for k = 1:numel(Zs) + Zs[k] = (D[k].bbox[5] + D[k].bbox[6]) / 2 + end + # But we need to check also if we have a vertical wall. In that case we set Zs[k] = 1e6 so that + # we can send these triangles to the foreground color of the CPT. + if (!isempty(D[1].comment) && contains(D[1].comment[1], "vwall")) + for k = 1:numel(Zs) + startswith(D[k].header, "W ") && (Zs[k] = 1e6) + end + end + return Zs +end diff --git a/src/utils.jl b/src/utils.jl index eb8adea37..bf62fa69d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1097,26 +1097,65 @@ end # ---------------------------------------------------------------------------- """ - n = facenorm(M::Matrix{<:Real}; normalize=true) + n = facenorm(M::Matrix{<:Real}; normalize=true, zfact=1.0) Calculates the normal vector of a polygon with vertices in `M`. - `normalize`: By default, it returns the unit vector. If `false`, it returns the non-normalized vector that represents the area of the polygon. +- `zfact`: If different from 1, the last dimension of `M` (normally, the Z variable) + is multiplied by `zfact`. + ### Returns A 3 elements vector with the components of the normal vector. """ -function facenorm(M::Matrix{<:Real}; normalize=true) +function facenorm(M::Matrix{<:Real}; zfact=1.0, normalize=true) # Must take care of the case when the first and last vertices are the same - last = (M[1,:] == M[end,:]) ? size(M,1)-1 : size(M,1) - c = cross(M[last,:], M[1,:]) # First edge - for k = 1:last-1 # Loop from first to end-1 - c += cross(M[k,:], M[k+1,:]) # Sum the rest of edges + p1 = M[1,:] + last = (p1 == M[end,:]) ? size(M,1)-1 : size(M,1) + if (zfact == 1.0) + c = cross(M[last,:], p1) # First edge + for k = 1:last-1 # Loop from first to end-1 + c += cross(M[k,:], M[k+1,:]) # Sum the rest of edges + end + else + p2 = M[last,:] + p1[end] *= zfact; p2[end] *= zfact + c = cross(p2, p1) # First edge + for k = 1:last-1 # Loop from first to end-1 + p1, p2 = M[k,:], M[k+1,:] + p1[end] *= zfact; p2[end] *= zfact + c += cross(p1, p2) # Sum the rest of edges + end end return normalize ? c / norm(c) : c # The cross product gives us 2 * A(rea) end +# --------------------------------------------------------------------------------------------------- +""" + A, B, C, D = eq_plane(azim, elev, dist) + +Calculate the equation of a plane (Eq: Ax + By + Cz + D = 0) given the azimuth, +elevation and distance to the plane. + +- `azim`: Azimuth in degrees (clockwise from North) +- `elev`: Elevation in degrees +- `dist`: Distance to the plane + +To compute the distance fom a point to that plane, do: +```julia + p = (0,0,0); + dist = abs(A * p[1] + B * p[2] + C * p[3] + D) +``` +""" +function eq_plane(azim, elev, dist) + nx, ny, nz = cosd(elev) * sind(azim), cosd(elev) * cosd(azim), sind(elev) + Px, Py, Pz = dist * nx, dist * ny, dist * nz + D = Px * nx + Py * ny + Pz * nz + return nx, ny, nz, D # Eq of plane: Ax + By + Cz + D = 0 +end + # --------------------------------------------------------------------------------------------------- include("makeDCWs.jl") @@ -1148,5 +1187,3 @@ end =# #GI.geometry[1].geoms[1].rings[1].vertices.data[1].coords.lat.val - -include("circfit.jl") \ No newline at end of file diff --git a/test/test_utils.jl b/test/test_utils.jl index 497585c6c..e17d8ea14 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -31,5 +31,7 @@ x2, y2, R2, err2 = circfit(x, y, taubin=true) return [x y], x1, y1, R1, err1, x2, y2, R2, err2 end - test_circfit(); + test_circfit(); + + A, B, C, D = GMT.eq_plane(0, 45, 10); end