Finish the grid2tri function that created 3D triangulated bodies from grids.

src/GMT.jl

@@ -144,13 +144,15 @@ export
 	getregion, gd2gmt, gmt2gd, gdalread, gdalshade, gdalwrite, gadm, xyzw2cube, coastlinesproj, graticules, orbits, orbits!,
 	plotgrid!, worldrectangular, worldrectgrid,
 
-	earthregions, gridit, magic, rescale, stackgrids, delrows!, setgrdminmax!, meshgrid, cart2pol, pol2cart, cart2sph, sph2cart,
+	earthregions, gridit, grid2tri, magic, rescale, stackgrids, delrows!, setgrdminmax!, meshgrid, cart2pol, pol2cart,
+	cart2sph, sph2cart,
 
 	arcellipse, arccircle, getband, getdriver, getlayer, getproj, getgeom, getgeotransform, gdaldrivers, toPROJ4, toWKT,
 	importPROJ4, importWKT, importEPSG, gdalinfo, gdalwarp, gdaldem, gdaltranslate, gdalgrid, gdalvectortranslate,
 	ogr2ogr, gdalrasterize, gdalbuildvrt, readgeom, readraster, setgeotransform!, setnodata!, setproj!, destroy,
 	delaunay, dither, buffer, centroid, intersection, intersects, polyunion, overlaps, fromWKT, fillnodata!, fillnodata,
 	concavehull, convexhull, difference, symdifference, distance, geodesicarea, geomarea, pointalongline, polygonize, simplify,
+
 	wkbUnknown, wkbPoint, wkbPointZ, wkbLineString, wkbPolygon, wkbPolygonZM, wkbMultiPoint, wkbMultiPointZ, wkbMultiLineString,
 	wkbMultiPolygon, wkbGeometryCollection, wkbPoint25D, wkbLineString25D, wkbPolygon25D, wkbMultiPoint25D,
 	wkbMultiLineString25D, wkbMultiPolygon25D, wkbGeometryCollection25D,

src/common_options.jl

@@ -1594,7 +1594,7 @@ function parse_helper(cmd::String, d::Dict, symbs::VMs, opt::String, sep='/')
 	(SHOW_KWARGS[1]) && return (print_kwarg_opts(symbs, "(Common option not yet expanded)"),"")
 	opt_val::String = ""
 	if ((val = find_in_dict(d, symbs, true)[1]) !== nothing)
-		opt_val = opt * arg2str(val, sep)::String
+		opt_val = opt * arg2str(val, sep)
 		cmd *= opt_val
 	end
 	return cmd, opt_val
@@ -1658,7 +1658,7 @@ function parse_common_opts(d::Dict, cmd::String, opts::VMs, first::Bool=true)
 		if (opt_p != "")
 			if (opt_p == " -pnone")  CURRENT_VIEW[1] = "";	cmd = cmd[1:end-7];	opt_p = ""
 			elseif (startswith(opt_p, " -pa") || startswith(opt_p, " -pd"))
-				CURRENT_VIEW[1] = " -p210/30";
+				CURRENT_VIEW[1] = " -p217.5/30";
 				cmd = replace(cmd, opt_p => "") * CURRENT_VIEW[1]::String		# auto, def, 3d
 			else
 				CURRENT_VIEW[1] = opt_p
@@ -1669,7 +1669,6 @@ function parse_common_opts(d::Dict, cmd::String, opts::VMs, first::Bool=true)
 			CURRENT_VIEW[1] = ""		# Ensure we start empty
 		end
 	end
-	#cmd   = GMTsyntax_opt(d, cmd)		# See if an hardcore GMT syntax string has been passed
 	((val = find_in_dict(d, [:pagecolor])[1]) !== nothing) && (cmd *= string(" --PS_PAGE_COLOR=", val)::String)
 	return cmd, o
 end

src/gdal_extensions.jl

@@ -37,14 +37,16 @@ end
 """
     buffer(geom, dist::Real, quadsegs::Integer=30; gdataset=false)
 
+Compute buffer of a geometry.
+
 ### Parameters
 * `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
 * `dist`: the buffer distance to be applied. Should be expressed into the
     same unit as the coordinates of the geometry.
 * `quadsegs`: the number of segments used to approximate a 90 degree (quadrant) of curvature.
-* `gdataset`: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix
 
-Compute buffer of geometry.
+### Keywords
+* `gdataset`: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix
 
 Builds a new geometry containing the buffer region around the geometry on which
 it is invoked. The buffer is a polygon containing the region within the buffer
@@ -73,6 +75,8 @@ end
 
 ### Parameters
 * `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of them), or a Matrix
+
+### Keywords
 * `gdataset`: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix
 
 Compute the geometry centroid.
@@ -107,6 +111,8 @@ end
 ### Parameters
 * `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
 * `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type 
+
+### Keywords
 * `gdataset`: Returns a GDAL IGeometry even when input is GMTdataset or Matrix
 
 Returns a new geometry representing the intersection of the geometries, or NULL
@@ -125,12 +131,14 @@ intersection(D1, D2; gdataset=false) = helper_geoms_run_fun(intersection, D1, D2
 """
     polyunion(geom1, geom2; gdataset=false)
 
+Computes a new geometry representing the union of the geometries.
+
 ### Parameters
 * `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
 * `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type 
-* `gdataset`: Returns a GDAL IGeometry even when input is GMTdataset or Matrix
 
-Computes a new geometry representing the union of the geometries.
+### Keywords
+* `gdataset`: Returns a GDAL IGeometry even when input is GMTdataset or Matrix
 
 ### Returns
 A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
@@ -145,6 +153,8 @@ polyunion(D1, D2; gdataset=false) = helper_geoms_run_fun(polyunion, D1, D2; gdat
 ### Parameters
 * `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
 * `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type 
+
+### Keywords
 * `gdataset`: Returns a GDAL IGeometry even when input is GMTdataset or Matrix
 
 Generates a new geometry which is the region of this geometry with the region of the other geometry removed.
@@ -162,6 +172,8 @@ difference(D1, D2; gdataset=false) = helper_geoms_run_fun(difference, D1, D2; gd
 ### Parameters
 * `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
 * `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type 
+
+### Keywords
 * `gdataset`: Returns a GDAL IGeometry even when input is GMTdataset or Matrix
 
 Generates a new geometry representing the symmetric difference of the geometries
@@ -267,6 +279,8 @@ Compute a simplified geometry.
 ### Parameters
 * `geom`: the geometry.
 * `tol`: the distance tolerance for the simplification.
+
+### Keywords
 * `gdataset`: Returns a GDAL IGeometry even when input is GMTdataset or Matrix
 """
 simplify(geom::AbstractGeometry, tol::Real) = IGeometry(OGR_G_Simplify(geom.ptr, tol))
@@ -486,6 +500,7 @@ end
 * `radius`: radius of the circle.
 * `start_angle`: angle to first point on arc (clockwise of X-positive) 
 * `end_angle`: angle to last point on arc (clockwise of X-positive) 
+
 ### Keywords
 * `z0`: center Z. Optional, if not provided the output is flat 2D
 * `inc`: the largest step in degrees along the arc. Default is 2 degree

src/gmtreadwrite.jl

@@ -706,3 +706,42 @@ function read_obj(fname)
 	set_dsBB!(DV)
 	return [DV, GMTdataset(data=F)]
 end
+
+# --------------------------------------------------------------------------------------------------------
+"""
+    write_stl(D::Vector{<:GMTdataset}, fname; binary=true)
+
+Write a STL file. 
+"""
+# Inspired in MeshIO stl.jl
+function write_stl(fname::AbstractString, D::Vector{<:GMTdataset}; binary::Bool=true)
+	name = fileparts(fname)[2]
+	fid = open(fname, write=true)
+
+	if (binary)
+		foreach(k -> write(fid, 0x00), 1:80)	# header (empty)
+		write(fid, UInt32(length(D)))			# number of triangles
+    	for k = 1:length(D)
+			n = facenorm(D[k].data)
+			foreach(j-> write(fid, Float32(n[j])), 1:3)
+			for t = 1:3
+				write(fid, Float32(D[k][t,1]), Float32(D[k][t,2]), Float32(D[k][t,3]))
+			end
+			write(fid, 0x0000)		# write 16bit empty byte count
+		end
+	else
+    	write(fid, "solid $name\n")
+		for k = 1:length(D)
+			n = facenorm(D[k].data)
+			@printf fid "facet normal %.12g %.12g %.12g\n" n[1] n[2] n[3]
+			write(fid,"\touter loop\n")
+			@printf fid "\t\tvertex  %.12g %.12g %.12g\n" D[k][1,1] D[k][1,2] D[k][1,3]
+			@printf fid "\t\tvertex  %.12g %.12g %.12g\n" D[k][2,1] D[k][2,2] D[k][2,3]
+			@printf fid "\t\tvertex  %.12g %.12g %.12g\n" D[k][3,1] D[k][3,2] D[k][3,3]
+			write(fid,"\tendloop\n")
+			write(fid,"endfacet\n")
+		end
+		write(fid,"endsolid $name\n")
+	end
+	close(fid)
+end

src/psxy.jl

@@ -174,6 +174,7 @@ function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D::
 	# Look for color request. Do it after error bars because they may add a column
 	len_cmd = length(cmd);	n_prev = N_args;
 	opt_Z, args, n, got_Zvars = add_opt(d, "", "Z", [:Z :level :levels], :data, Any[arg1, arg2], (outline="_o", nofill="_f"))
+	(opt_Z == " -Z" && n == 0) && error("The 'level' option (Z) must be set a single value, a file name or a vector os reals.")
 	if (isa(arg1, Vector{<:GMTdataset}) && ((ind_att = findfirst('=', opt_Z)) !== nothing))
 		# Here we deal with the case where the -Z option refers to a particular attribute.
 		# Allow to use "Z=(data="att=XXXX", nofill=true)" when the Di are polygons.

src/triangulate.jl

@@ -174,20 +174,155 @@ end
 trisurf!(in::Union{Matrix, GDtype}; gdal=true, kw...) = trisurf(in; gdal=gdal, first=false, kw...)
 
 # ---------------------------------------------------------------------------------------------------
-function grd2FV(G)
-	isa(G, String) && (G = gmtread(G))
-	#V = grd2xyz(G, s=true)
-	#V.geom = wkbPointZ
-	#F = Int.(triangulate(V, T=true).data) .+ 1
-	#return [V, GMTdataset(data=F)]
-	V = grd2xyz(G, s=true)
-	T = triangulate(V, S=true, Z=true)
-	C = gmtspatial(T, Q=true, o="0,1")
-	#C = centroid(T)
-	B = concavehull(V.data, 0.001)
-	ind = (C in B) .== 1
-	T = T[ind]
-	T[1].ds_bbox = T[1].bbox			# Because we may have deleted first T
-	set_dsBB!(T, false)
-	return T
+# D = grid2tri("cam_slab2_dep_02.24.18.grd", "cam_slab2_thk_02.24.18.grd")
+# plot3(D,  zsize=3, G="+z", L=true, proj=:merc, show=1, pen=0)
+# plot3(D,  zsize=3,  proj=:merc,  G="+z", C="aa.cpt", show=1, Z=linspace(-7,-400,length(D)), Vd=1)
+"""
+    grid2tri(G, G2=nothing; bottom=false, downsample=0, level=false, ratio=0.01, thickness=0.0, wall=false)
+
+Triangulates the surface defined by the grid `G`, and optionally the bottom surface `G2`, plus the
+vertical wall between them, or between `G` and constant level or a constant thickness. Optionally
+computes only the vertical wall or the full closed bodie (that is, including the bottom surface).
+
+The output of this function can be used in ``plot3d`` to create 3D views of volume layer.
+
+NOTE: The `G` grid should have a _out-skirt_ of NaNs, otherwise just use ``grdview`` with the ``N`` option.
+
+### 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
+  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.
+
+- `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.
+
+- `downsample`: If the grid is of too high resolution, files here get big and operations slow down with this
+  and later figures may not benefit much. In those cases it is a good idea to downsample the grid. The 
+  `downsample` option accepts an integer reduction factor. `downsample=2` will shrink the grid by a factor two
+  in each dimention, `downsample=3` will shrink it by a factor three etc.
+
+- `ratio`: A slightly tricky parameter that determines how close the computed concave hull is to the true
+  concave hull. A value smaller to 0.005 seems to do it but we normally don't want that close because the
+  vertical wall obtained from this will be too jagged. The default value of 0.01 seems to work well to get
+  a smoother concave hull but one that still fits the objective of getting a nice vertical wall. May need
+  tweaking for specific cases.
+
+- `thickness`: A scalar representing the layer thickness in the same units as those of the input grid.
+  NOTE: this option is ignored when two grids are passed in input.
+
+- `wall`: If true, only the vertical wall  between `G` and `G2`, or `G` + `thickness` is computed and returned.
+
+### 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, wall=false)
+	(!isa(G2, GMTgrid) && !isa(G2, String) && thickness <= 0.0) &&
+		error("Must provide a bottom grid or a thickness for this layer.")
+
+	(wall != 0) && (bottom = false)
+	Dbnd_t, Dpts = gridhull(G; downsample=downsample, ratio=ratio)	# Compute the top concave hull
+	if (wall != 0)										# If we only want the vertical wall no need to compute the others
+		Dt_t = triangulate(Dpts, S="+za", 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[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
+		set_dsBB!(Dt_t, false)
+	end
+
+	if (isa(G2, GMTgrid) || isa(G2, String))
+		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 = Dt_b[ind]							# Delete the triangles outside the concave hull (reuse the same 'ind')
+			Dt_b[1].ds_bbox = Dt_b[1].bbox				# Because we may have deleted first Dt
+			Dt_b[1].geom = wkbPolygonZM
+			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))
+		(bottom == 1) && append!(Dt_b, Dwall)			# If including bottom too, start with it and add the wall.
+	else
+		Dwall = vwall(Dbnd_t, thickness, level != 0)
 	end
+
+	(wall == 0) && append!(Dwall, Dt_t)
+	return Dwall
+end
+
+# ---------------------------------------------------------------------------------------------------
+"""
+    B, V = gridhull(G; downsample::UInt=0, ratio=0.01) -> GMTdataset, Matrix
+
+- `G`: The input grid. It can be either a GMTgrid or a grid file name.
+
+### Keywords
+- `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
+	(downsample > 1) && (_G = grdsample(_G, I="$(div(size(_G.z,2),2)+1)+n/$(div(size(_G.z,1),2)+1)+n", V="q"))
+	V = grd2xyz(_G, s=true)					# Convert to x,y,z while dropping the NaNs
+	B = concavehull(V.data, ratio, false)	# Compute the ~concave hull when excluding the outer NaNs (ignores holes)
+	B.proj4, B.wkt, B.epsg = _G.proj4, _G.wkt, _G.epsg
+	return B, V
+end
+
+# ---------------------------------------------------------------------------------------------------
+"""
+    D = vwall(Bt, thk, level::Bool=false) -> Vector{GMTdataset}
+
+Compute the vertical wall between grid's concave hull `Bt` with a fixed or variable height `thk`.
+
+- `Bt`: A Mx3 matrix or a GMTdataset with the concave hull of the top surface given in a clock-wise order
+  (The order returned by GDAL's `concavehull` function).
+
+- `thk`: A constant or a vector with the thickness of the wall.
+
+- `level`: If ``true`` `thk` is interpreted as the level of the bottom surface instead of a constant thickness.
+"""
+function vwall(Bt::Union{Matrix{<:Real}, GMTdataset}, thk::Union{<:Real, AbstractVector{<:Real}}, level::Bool=false)
+	Bb = copy(Bt)						# Always return a Matrix, not a DS
+	if isa(thk, Real)
+		if (level)  view(Bb, :, 3) .=  convert(eltype(Bt), thk)
+		else        view(Bb, :, 3) .-= convert(eltype(Bt), thk)
+		end
+	else            view(Bb, :, 3) .-= thk
+	end
+	vwall(Bt, Bb)
+end
+
+# ---------------------------------------------------------------------------------------------------
+"""
+    D = vwall(Bt::Union{Matrix, GMTdataset}, Bb::Union{Matrix, GMTdataset}) -> Vector{GMTdataset}
+
+Compute the vertical wall between two grid's concave hull `Bt` and `Bb`.
+"""
+function vwall(Bt::Union{Matrix{<:Real}, GMTdataset}, Bb::Union{Matrix{<:Real}, GMTdataset})
+	(size(Bt) != size(Bb)) && error("Input matrices must be the same size")
+	n_sideT = size(Bb,1) - 1
+	Twall = Vector{GMTdataset}(undef, 2 * n_sideT)
+	for k = 1:n_sideT
+		kk = 2 * k -1
+		Twall[kk]   = GMTdataset(data=[Bt[k,1]   Bt[k,2]   Bt[k,3];   Bt[k+1,1] Bt[k+1,2] Bt[k+1,3]; Bb[k,1] Bb[k,2] Bb[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]])
+		Twall[kk].header = Twall[kk+1].header = "-Z$(Bt[k,3]+1e6)"		# Can be used, e.g., with the foreground color of a CPT 
+	end
+	set_dsBB!(Twall)
+	isa(Bt, GMTdataset) && (Twall[1].proj4 = Bt.proj4)
+	Twall[1].geom = wkbPolygonZM
+	return Twall
+end