Merge pull request #1369 from GenericMappingTools/read-cubes-with-V

Let gdalread also fetch the v coordinates for cubes.

src/GMT.jl

@@ -166,7 +166,7 @@ export
 
 	ablines, ablines!, density, density!, boxplot, boxplot!, cornerplot, cornerplot!, cubeplot, ecdfplot, ecdfplot!,
 	fill_between, fill_between!, marginalhist, marginalhist!, parallelplot, parallelplot!, plotlinefit, plotlinefit!,
-	qqplot, qqplot!, qqnorm, qqnorm!, seismicity, sealand, terramar, violin, violin!, viz, weather, windbarbs, whereami,
+	qqplot, qqplot!, qqnorm, qqnorm!, seismicity, sealand, squeeze, terramar, violin, violin!, viz, weather, windbarbs, whereami,
 
 	VSdisp, info, kmeans, pca, mosaic, quadbounds, quadkey, geocoder,
 

src/gdal_utils.jl

@@ -109,9 +109,18 @@ function gd2gmt(_dataset; band::Int=0, bands=Vector{Int}(), sds::Int=0, pad::Int
 	(prj != "" && !startswith(prj, "+proj")) && (prj = toPROJ4(importWKT(prj)))
 	(prj == "") && (prj = seek_wkt_in_gdalinfo(gdalinfo(dataset)))
 	is_tp = (layout == "")				# if == "" array is rowmajor and hence transposed
+	vname, vvalues = "", Float64[]
+	if (n_dsbands > 1 && ((meta = Gdal.getmetadata(dataset)) != String[]))	# For cubes, try to find the v coordinates
+		if ((ind = findfirst(startswith.(meta, "NETCDF_DIM_EXTRA="))) !== nothing)
+			vname = meta[ind][19:end-1]		# The variable name is wrapped in {}
+			if ((ind = findfirst(startswith.(meta, "NETCDF_DIM_" * vname * "_VALUES="))) !== nothing)
+				vvalues = parse.(Float64, split(meta[ind][length(vname)+21:end-1], ","))	# for ex: NETCDF_DIM_Depth_VALUES={
+			end
+		end
+	end
 	if (is_grid)
 		#(eltype(mat) == Float64) && (mat = Float32.(mat))
-		O = mat2grid(mat; hdr=hdr, proj4=prj, names=desc, is_transposed=is_tp)
+		O = mat2grid(mat; hdr=hdr, v=vvalues, v_unit=vname, proj4=prj, names=desc, is_transposed=is_tp)
 		O.layout = (layout == "") ? "TRB" : layout
 		isa(mat, Matrix{<:Complex}) && (append!(O.range, [z_im_min, z_im_max]))	# Stick the imaginary part limits in the range
 	else

src/spectrum1d.jl

@@ -21,7 +21,7 @@ Parameters
 - **L** | **leave_trend** :: [Type => Str | []]     ``Arg = [h|m]``
 
     Leave trend alone. By default, a linear trend will be removed prior to the transform.
-- **N** | **time_col** :: [Type => Int]      ``Arg = t_col``
+- **N** | **name** :: [Type => Int]      ``Arg = t_col``
 
     Indicates which
 - **T** :: [Type => Bool]
@@ -52,7 +52,7 @@ function spectrum1d_helper(cmd0::String, arg1; kwargs...)
 
 	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
 	cmd = parse_common_opts(d, "", [:V_params :b :d :e :g :h :i :yx])[1]
-	cmd = parse_these_opts(cmd, d, [[:D :sample_dist], [:L :leave_trend], [:N :time_col], [:S :size], [:T :multifiles], [:W :wavelength]])
+	cmd = parse_these_opts(cmd, d, [[:D :sample_dist], [:L :leave_trend], [:N :name], [:S :size], [:T :multifiles], [:W :wavelength]])
 	opt_C = add_opt(d, "", "C", [:C :components :outputs :response_fun],
 	                (xpower="_x", ypower="_y", cpower="_c", npower="_n", phase="_p", admit="_a", gain="_g", coh="_o"))
     if (!contains(cmd, " -T"))                  # Return vars, otherwise it saves them in disk files

src/sphdistance.jl

@@ -14,7 +14,7 @@ Parameters
 - **D** | **duplicates** :: [Type => Bool]
 
     Delete any duplicate points [Default assumes there are no duplicates].
-- **E** | **what_quantity** :: [Type => Str]   ``Arg = d|n|z[dist]``
+- **E** | **quantity** :: [Type => Str]   ``Arg = d|n|z[dist]``
 
     Specify the quantity that should be assigned to the grid nodes.
 - **G** | **save** | **outgrid** | **outfile** :: [Type => Str]
@@ -49,7 +49,7 @@ function sphdistance(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)
 	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
 
 	cmd, = parse_common_opts(d, "", [:G :RIr :V_params :b :d :e :h :i :yx])
-	cmd  = parse_these_opts(cmd, d, [[:C :save_mem], [:D :duplicates], [:E :what_quantity], [:L :dist_unit]])
+	cmd  = parse_these_opts(cmd, d, [[:C :save_mem], [:D :duplicates], [:E :quantity], [:L :dist_unit]])
 	cmd, arg1, arg2 = parse_QN_sphdst(d, [[:Q :voronoi], [:N :nodes]], cmd, arg1, arg2)
 
 	common_grd(d, "sphdistance " * cmd, arg1, arg2)