Rework the 3D bodies code to work with the new GMTfv type. (#1559)

* Remove temporary wrappers.

* Identify the new GMTfv type

* When buffergeo is given single point, don't try to do distance to line.

* Remove the now obsolete isFV function

* DEclare and export the new GMTfv type.

* Chande the dataset_init_FV to deal with the new GMTfv type

* Add fv2fv function and methods to create instances of GMTfv

* Be aware of GMTfv types

* Return the solids as GMTfv types.

* Rework the 3D bodies code to work with the new GMTfv type.

Author: joa-quim

src/GMT.jl

@@ -112,7 +112,7 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optle
 end
 
 export
-	GMTgrid, GMTimage, GMTdataset, GMTcpt, GItype, GDtype, GMTver, FMT, TMPDIR_USR, gmt, libgdal,
+	GMTgrid, GMTimage, GMTdataset, GMTfv, GMTcpt, GItype, GDtype, GMTver, FMT, TMPDIR_USR, gmt, libgdal,
 	arrows, arrows!, bar, bar!, bar3, bar3!, band, band!, bubblechart, bubblechart!, feather, feather!, hband, hband!,
 	hlines, hlines!, lines, lines!, legend, legend!, quiver, quiver!, radar, radar!, stairs, stairs!, stem, stem!,vlines,
 	vlines!, vband, vband!, hspan, hspan!, vspan, vspan!,
@@ -186,8 +186,7 @@ export
 
 	cube, dodecahedron, icosahedron, sphere, octahedron, tetrahedron, replicant,
 
-	df2ds, ODE2ds,
-	sprintf,
+	df2ds, ds2df, ODE2ds,
 	@?, @dir
 
 include("common_docs.jl")

src/common_options.jl

@@ -3643,6 +3643,14 @@ end
 # ---------------------------------------------------------------------------------------------------
 function read_data(d::Dict, fname::String, cmd::String, arg, opt_R::String="", is3D::Bool=false, get_info::Bool=false)
 	# Use 'get_info=true' to force reading the file when fname != ""
+	
+	if (isa(arg, GMTfv))			# A quick and dirty way to parse the GMTfv type
+		is_geo = isgeog(arg.proj4)			# MUST make isgeog sniff in GMTfv types
+		_wesn = round_wesn(arg.bbox, is_geo)		# Add a pad if not-tight
+		opt_R = @sprintf(" -R%.12g/%.12g/%.12g/%.12g/%.12g/%.12g", _wesn[1], _wesn[2], _wesn[3], _wesn[4], _wesn[5], _wesn[6])
+		return cmd * opt_R, arg, opt_R, _wesn, ""		# A GMTfv is already read
+	end
+
 	cmd::String, opt_i::String  = parse_i(d, cmd)	# If data is to be read with some column order
 	cmd, opt_bi::String = parse_bi(d, cmd)			# If data is to be read as binary
 	cmd, opt_di::String = parse_di(d, cmd)			# If data missing data other than NaN
@@ -3660,7 +3668,7 @@ function read_data(d::Dict, fname::String, cmd::String, arg, opt_R::String="", i
 		if (opt_yx != "" && !isempty(arg.colnames))  arg.colnames[2], arg.colnames[1] = arg.colnames[1], arg.colnames[2]   end
 		# Try guess if ascii file has time columns and if yes leave trace of it in GMTdadaset metadata.
 		(opt_bi == "") && file_has_time!(fname, arg)	# If fname is a .gmt file this does not make much sense.
-		# Remove the these options from cmd. Their job is done
+		# Remove these options from cmd. Their job is done
 		if (opt_i != "")  cmd = replace(cmd, opt_i => "");	opt_i = ""  end
 		if (opt_h != "")  cmd = replace(cmd, opt_h => "");	opt_h = ""  end
 	else							# No need to find -R so let the GMT module read the file
@@ -3714,9 +3722,9 @@ function _read_data(d::Dict, cmd::String, arg, opt_R::String="", is3D::Bool=fals
 
 	have_info = false
 	no_R = (opt_R == "" || opt_R[1] == '/' || opt_R == " -Rtight")
-	prj::String = (isGMTdataset(arg)) ? getproj(arg, proj4=true) : ""
+	prj::String = (isa(arg, GDtype)) ? getproj(arg, proj4=true) : ""
 	is_geo = isgeog(prj)
-	ds_bbox = isGMTdataset(arg) ? (isa(arg, GMTdataset) ? arg.ds_bbox : arg[1].ds_bbox) : Float64[]
+	ds_bbox = isa(arg, GDtype) ? (isa(arg, GMTdataset) ? arg.ds_bbox : arg[1].ds_bbox) : Float64[]
 
 	if (no_R && !isempty(ds_bbox))			# If arg is a GMTdataset of polygons, use the ds_bbox directly
 		geom = isa(arg, Vector) ? arg[1].geom : arg.geom
@@ -3943,17 +3951,6 @@ function isgeog(in)::Bool
 	(prj != "" && (contains(prj, "=lon") || contains(prj, "=lat")))
 end
 
-# ---------------------------------------------------------------------------------------------------
-"""
-    isFV(D)::Bool
-
-Check if D is a Face-Vertices ensemble (a 2 elements vector of GMTdataset).
-"""
-function isFV(D)::Bool
-	(isa(D, Vector{<:GMTdataset}) && length(D) > 1) &&
-		(D[1].geom == wkbPoint || D[1].geom == wkbPointZ) && eltype(D[2].data) <: Integer
-end
-
 # ---------------------------------------------------------------------------------------------------
 """
     is_gridtri(D)::Bool

src/gmt_main.jl

@@ -722,6 +722,8 @@ function GMTJL_Set_Object(API::Ptr{Nothing}, X::GMT_RESOURCE, ptr, pad)::GMT_RES
 			end
 		elseif (isa(ptr, Vector{<:GMTdataset}))
 			X.object = dataset_init(API, ptr, X.direction)
+		elseif (isa(ptr, GMTfv))
+			X.object = dataset_init_FV(API, ptr)
 		else
 			if (X.direction == GMT_OUT)		# Here we accept ptr === nothing
 				X.object = convert(Ptr{GMT_DATASET}, GMT_Create_Data(API, GMT_IS_DATASET, GMT_IS_PLP, GMT_IS_OUTPUT, NULL, NULL, NULL, 0, 0, NULL))
@@ -1017,8 +1019,6 @@ function dataset_init(API::Ptr{Nothing}, Darr::Vector{<:GMTdataset}, direction::
 
 	(Darr == C_NULL || length(Darr) == 0) && error("Input is empty where it can't be.")
 
-	isFV(Darr) && return dataset_init_FV(API, Darr)			# Special case for Face-Vertices obj.
-
 	# We come here if we did not receive a matrix
 	dim = [1, 0, 0, 0]
 	dim[GMT_SEG+1] = length(Darr)					# Number of segments
@@ -1091,33 +1091,31 @@ end
 
 # ---------------------------------------------------------------------------------------------------
 function dataset_init_FV(API::Ptr{Nothing}, FV)::Ptr{GMT_MATRIX}
-	V::GMTdataset{Float64,2}, F::GMTdataset{Int,2} = FV[1], FV[2]
-	n_segs = size(F.data, 1)				# Number of segments or faces (polygons)
-	n_rows = size(F.data, 2)				# Number of rows (vertexes of the polygon)
-	n_cols = size(V.data, 2)				# Number of columns (2 for x,y; 3 for x,y,z)
-	dim = [1, n_segs, n_rows, n_cols]		# [1, GMT_SEG+1, GMT_ROW+1, GMT_COL+1]
+	n_segs = size(FV.faces[1], 1)			# Number of segments or faces (polygons)
+	n_rows = size(FV.faces[1], 2)			# Number of rows (vertices of the polygon)
+	dim = [1, n_segs, n_rows, 3]			# [1, GMT_SEG+1, GMT_ROW+1, GMT_COL+1]
 
 	pdim = pointer(dim)
 	D = convert(Ptr{GMT_DATASET}, GMT_Create_Data(API, GMT_IS_DATASET, GMT_IS_PLP, GMT_NO_STRINGS, pdim, NULL, NULL, 0, 0, NULL))
 	DS::GMT_DATASET = unsafe_load(D)
-	DT = unsafe_load(unsafe_load(DS.table))			# GMT_DATATABLE
+	DT = unsafe_load(unsafe_load(DS.table))		# GMT_DATATABLE
 
 	n_records = 0
-	tmp = zeros(n_rows, n_cols)
+	tmp = zeros(n_rows, 3)
 
-	for seg = 1:n_segs 								# Each row in F (a face) is a new data segment (a polygon)
+	for seg = 1:n_segs 							# Each row in a face is a new data segment (a polygon)
 		DSv = convert(Ptr{Nothing}, unsafe_load(DT.segment, seg))		# DT.segment = Ptr{Ptr{GMT_DATASEGMENT}}
-		S = GMT_Alloc_Segment(API, GMT_NO_STRINGS, n_rows, n_cols, "", DSv) # Ptr{GMT_DATASEGMENT}
-		Sb = unsafe_load(S)							# GMT_DATASEGMENT;		Sb.data -> Ptr{Ptr{Float64}}
-		
-		for c = 1:n_cols, r = 1:n_rows
-			tmp[r,c] = V.data[F.data[seg, r], c]
+		S = GMT_Alloc_Segment(API, GMT_NO_STRINGS, n_rows, 3, "", DSv) # Ptr{GMT_DATASEGMENT}
+		Sb = unsafe_load(S)						# GMT_DATASEGMENT;		Sb.data -> Ptr{Ptr{Float64}}
+	
+		for c = 1:3, r = 1:n_rows
+			tmp[r,c] = FV.verts[FV.faces[1][seg, r], c]
 		end
-		for col = 1:n_cols							# Copy the data columns
+		for col = 1:3							# Copy the data columns
 			unsafe_copyto!(unsafe_load(Sb.data, col), pointer(tmp[:,col]), n_rows)
 		end
 
-		n_records += n_rows							# Must manually keep track of totals
+		n_records += n_rows						# Must manually keep track of totals
 		DS.type_ = GMT_READ_DATA
 		unsafe_store!(S, Sb)
 		unsafe_store!(DT.segment, S, seg)
@@ -1579,6 +1577,11 @@ end
 Base.:show(io::IO, ::MIME"text/plain", D::GMTdataset) = show(D)
 Base.:display(D::GMTdataset) = show(D)		# Otherwise the default prints nothing when text only (data == [])
 
+# ---------------------------------------------------------------------------------------------------
+function Base.:show(io::IO, ::MIME"text/plain", FV::GMTfv)
+	println("Face-Vertices with a total of ", sum(size.(FV.faces,1)), " faces and $(size(FV.verts,1)) vertices")
+end
+
 # ---------------------------------------------------------------------------------------------------
 function Base.show(io::IO, C::GMTcpt)
 	isempty(C) && return

src/gmt_types.jl

@@ -292,6 +292,17 @@ GMTdataset(data::Array{Float32,2}, text::String) =
 GMTdataset(data::Array{Float32,2}) =
 	GMTdataset(data, Float64[], Float64[], DictSvS(), String[], String[], "", String[], "", "", 0, 0)
 
+Base.@kwdef struct GMTfv{T<:AbstractFloat} <: AbstractMatrix{T}
+	verts::AbstractMatrix{T}=Matrix{Float64}(undef,0,0)
+	faces::Vector{<:AbstractMatrix{Int}}=Vector{<:AbstractMatrix{Int}}(undef,0)
+	faces_view::Vector{Matrix{Int}}=Vector{Matrix{Int}}(undef,0)
+	bbox::Vector{Float64}=zeros(6)
+	proj4::String=""
+	wkt::String=""
+	epsg::Int=0
+end
+Base.size(FV::GMTfv) = sum(size.(FV.faces))
+
 #=
 Base.@kwdef struct GMTtypes
 	stored::String = ""

src/gmtreadwrite.jl

@@ -679,7 +679,7 @@ end
 
 # --------------------------------------------------------------------------------------------------------
 """
-    VF = read_obj(fname)
+    VF = read_obj(fname)::GMTfv
 
 Read a Wavefront .obj file and return the result in a FaceVertices object.
 """
@@ -725,9 +725,7 @@ function read_obj(fname)
 		end
 	end
 	close(fid)
-	DV = GMTdataset(data=V, geom=wkbPointZ)
-	set_dsBB!(DV)
-	return [DV, GMTdataset(data=F)]
+	fv2fv(F, V)
 end
 
 # --------------------------------------------------------------------------------------------------------

src/imshow.jl

@@ -65,14 +65,14 @@ function imshow(arg1, x::AbstractVector{Float64}=Float64[], y::AbstractVector{Fl
 	elseif (isa(arg1, Array{UInt8}) || isa(arg1, Array{UInt16,3}))
 		Gi = mat2img(arg1; kw...)
 		call_img = true
-	elseif (isGMTdataset(arg1) || (isa(arg1, Matrix{<:Real}) && size(arg1,2) <= 4) || (isa(arg1, Gdal.AbstractDataset) || isa(arg1, Gdal.AbstractGeometry)))
+	elseif (isa(arg1, GDtype) || isa(arg1, GMTfv) || (isa(arg1, Matrix{<:Real}) && size(arg1,2) <= 4) || (isa(arg1, Gdal.AbstractDataset) || isa(arg1, Gdal.AbstractGeometry)))
 		(isa(arg1, Gdal.AbstractDataset) || isa(arg1, Gdal.AbstractGeometry)) && (arg1 = gd2gmt(arg1))
 		isa(arg1, Matrix{<:Real}) && (arg1 = mat2ds(arg1))
-		ginfo = isa(arg1, GMTdataset) ? arg1.bbox : arg1[1].ds_bbox
+		ginfo = (isa(arg1, GMTdataset) || isa(arg1, GMTfv)) ? arg1.bbox : arg1[1].ds_bbox
 		CTRL.limits[1:4] = ginfo[1:4];		CTRL.limits[7:10] = ginfo[1:4]
 		call_plot3 = ((isa(arg1, GMTdataset) && arg1.geom == Gdal.wkbPolygonZM) ||
 		              (isa(arg1, Vector{<:GMTdataset}) && arg1[1].geom == Gdal.wkbPolygonZM) ||
-					  isFV(arg1)) ? true : false		# Should evolve into a fun that detects the several plot3d cases.
+					  isa(arg1, GMTfv)) ? true : false		# Should evolve into a fun that detects the several plot3d cases.
 		!call_plot3 && (call_plot3 = isplot3(kw))
 		return (call_plot3) ? plot3d(arg1; show=see, kw...) : plot(arg1; show=see, kw...)
 	elseif (isa(arg1, GMTcpt))

src/libgmt.jl

@@ -418,14 +418,6 @@ function gmt_get_rgb_from_z(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, value, rgb)
 	ccall((:gmt_get_rgb_from_z, libgmt), Cint, (Cstring, Ptr{Cvoid}, Cdouble, Ptr{Cdouble}), GMT_Get_Ctrl(API), P, value, rgb)
 end
 
-function gmt_get_index(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, value)
-	ccall((:gmt_get_index, libgmt), Cint, (Cstring, Ptr{Cvoid}, Ptr{Cdouble}), GMT_Get_Ctrl(API), P, value)
-end
-
-function gmt_get_rgb_lookup(API::Ptr{Cvoid}, P::Ptr{GMT_PALETTE}, ind, value, rgb)
-	ccall((:gmt_get_rgb_lookup, libgmt), Cint, (Cstring, Ptr{Cvoid}, Cint, Cdouble, Ptr{Cdouble}), GMT_Get_Ctrl(API), P, ind, value, rgb)
-end
-
 function gmt_illuminate(API::Ptr{Cvoid}, intensity, rgb::Vector{Float64})
 	ccall((:gmt_illuminate, libgmt), Cvoid, (Cstring, Cdouble, Ptr{Cdouble}), GMT_Get_Ctrl(API), intensity, rgb)
 end

src/plot.jl

@@ -1811,6 +1811,7 @@ const psevents = events            # Alias
 cat_1_arg(arg::GMTdataset, toDS::Bool=false) = return arg			# Miserable attempts to force type stability
 cat_1_arg(arg::Vector{<:GMTdataset}, toDS::Bool=false) = return arg
 function cat_1_arg(arg, toDS::Bool=false)
+	isa(arg, GMTfv) && return arg				# A FV type. Nothing to here
 	# Add a first column with 1:n to all args that are not GMTdatasets
 	if (isa(arg, Vector) || typeof(arg) <: AbstractRange)
 		if isa(arg, Vector{<:Vector{<:Real}})

src/proj_utils.jl

@@ -323,6 +323,10 @@ function buffergeo(line::Matrix{<:Real}; width=0, unit=:m, np=120, flatstart=fal
 	#(line[1, 1:2] == line[end, 1:2]) && (flatstart = flatend = false)		# Polygons can't have start/end flat edges
 	width *= unit_factor(unit)
 	n_seg = length(azim)
+	
+	# --------------------- When 'line' is actually a single point -------------------------------------------- 
+	n_seg == 0 && return mat2ds(circgeo(line, radius=width, np=np, proj=proj, epsg=epsg), proj=proj, epsg=epsg)
+
 	D, _D = GMTdataset(), GMTdataset()
 	for n = 1:n_seg
 		seg = [geod(line[n,:], azim[n], 0:width/4:dist[n], proj=proj, epsg=epsg)[1]; line[n+1:n+1,:]]