Merge pull request #1550 from GenericMappingTools/gridtri-more

Rework the grid2tri code to drop invisble triangles and reorder them so furthes triangles are plotted first.

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)))

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.

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...)
 	# ...

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

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

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
 

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