Move some functions around and replace some calls to module names by gmt(...)

src/GMT.jl

@@ -282,6 +282,7 @@ include("grdview.jl")
 include("grdvolume.jl")
 include("greenspline.jl")
 include("gridit.jl")
+include("histo_funs.jl")
 include("img_funs.jl")
 include("imgtiles.jl")
 include("imshow.jl")
@@ -291,6 +292,7 @@ include("loxodromics.jl")
 include("makecpt.jl")
 include("mapproject.jl")
 include("maregrams.jl")
+include("marker_name.jl")
 include("movie.jl")
 include("nearneighbor.jl")
 include("pastplates.jl")
@@ -345,6 +347,7 @@ include("extras/anaglyph.jl")
 include("extras/hampel_outliers.jl")
 include("extras/isoutlier.jl")
 include("extras/lowess.jl")
+include("extras/mapsize2region.jl")
 include("extras/seismicity.jl")
 include("extras/okada.jl")
 include("extras/weather.jl")

src/beziers.jl

@@ -62,15 +62,17 @@ function bezier(D::GMTdataset; t=nothing, np::Int=0, pure=false, firstcurve=true
 	is_geo = isgeog(D)
 	if (is_geo)
 		mat::Matrix{Float64} = (size(D.data, 2) == 2) ? [D.data fill(0, size(D.data, 1))] : D.data
-		mat = mapproject(mat, E=true).data						# Convert to ECEF
+		#mat = mapproject(mat, E=true).data						# Convert to ECEF
+		mat = gmt("mapproject -E", mat).data					# Convert to ECEF
 	else
 		mat = D.data
 	end
 	(pure == 1 && size(mat, 1) > 4) && (pure = false;
 		@warn("pure option not allowed for more than 4 control points. Reverting to 'impure'"))
 
 	out::Matrix{Float64} = bezier(mat; t=t, np=np, pure=pure, firstcurve=firstcurve)
-	is_geo && (out = mapproject(out, E=true, I=true).data)		# Convert back to geographic
+	#is_geo && (out = mapproject(out, E=true, I=true).data)		# Convert back to geographic
+	is_geo && (out = gmt("mapproject -E -I", out).data)			# Convert back to geographic
 	return mat2ds(out, D)::GMTdataset{Float64, 2}
 end
 

src/common_options.jl

@@ -4627,7 +4627,8 @@ function regiongeog(GI::GItype)::Tuple
 end
 function regiongeog(fname::String)::Tuple
 	((prj = getproj(fname, wkt=true)) == "") && (@warn("Input grid/image has no projection info"); return ())
-	info = grdinfo(fname, C=true);		# It should also report the
+	#info = grdinfo(fname, C=true);		# It should also report the
+	info = gmt("grdinfo -C " * fname);
 	c = xy2lonlat([info.data[1] info.data[3]; info.data[2] info.data[4]]; s_srs=prj)
 	tuple(c...)
 end
@@ -4961,7 +4962,8 @@ Check if `cmd0` is a remote grid or a OceanColor one and return the default CPT
 """
 function check_remote_cpt(cmd0::String)
 	(cmd0 == "") && return ""
-	cpt_path = joinpath(dirname(pathof(GMT)), "..", "share", "cpt")
+	#cpt_path = joinpath(dirname(pathof(GMT)), "..", "share", "cpt")
+	cpt_path = joinpath(dirname(pathof(getfield(Main, nameof(@__MODULE__)))), "..", "share", "cpt")
 	if     (occursin("SST.sst", cmd0))     return cpt_path * "/sst_oc.cpt"
 	elseif (occursin("CHL.chlor_a", cmd0)) return cpt_path * "/chlor_oc.cpt"
 	end

src/extras/anaglyph.jl

@@ -111,13 +111,19 @@ end
 # ---------------------------------------------------------------------------------------------------
 function anaglyph_3d(G::GMTgrid; zsize=4, azim=190, dazim=2, cmap="gray")
 
+	# We do this more convoluted call to grdview to avoid calling it directly (unknown in future GMT_base module)
 	p_left, p_right = "$(azim)/30/0", "$(azim - dazim)/30/0"
-	pato = TMPDIR_USR[1] * "/" * "GMTjl_" * TMPDIR_USR[2] * TMPDIR_USR[3]
-	grdview(G, J="Q20c", JZ="$zsize", B=:none, C=string(cmap), Q=:i, p=p_left, I=true, figname=pato * ".jpg")
-	Il = gdalread(pato * ".jpg", "-b 1", layout="TCBa")
+	pato = TMPDIR_USR[1] * "/" * "GMTjl_" * TMPDIR_USR[2] * TMPDIR_USR[3] * ".jpg"
+	opt_R = @sprintf("%.15g/%.15g/%.15g/%.15g/%.15g/%.15g", G.range[1], G.range[2], G.range[3], G.range[4], G.range[5], G.range[6])
+	#grdview(G, J="Q20c", JZ="$zsize", B=:none, C=string(cmap), Q=:i, p=p_left, I=true, figname=pato)
+	cmd = ["grdview -R$opt_R -JQ20c -JZ$zsize -C$cmap -Qi -p$p_left -I+a-45+nt1"]
+	prep_and_call_finish_PS_module(Dict{Symbol, Any}(:figname => pato), cmd, "", true, false, true, G)
+	Il = gdalread(pato, "-b 1", layout="TCBa")
 	nr_l, nc_l = size(Il,1), size(Il,2)
-	grdview(G, J="Q20c", JZ="$zsize", B=:none, C=string(cmap), Q=:i, p=p_right, I=true, figname=pato * ".jpg")
-	Ir = gdalread(pato * ".jpg", layout="TCBa")
+	#grdview(G, J="Q20c", JZ="$zsize", B=:none, C=string(cmap), Q=:i, p=p_right, I=true, figname=pato)
+	cmd = ["grdview -R$opt_R -JQ20c -JZ$zsize -C$cmap -Qi -p$p_right -I+a-45+nt1"]
+	prep_and_call_finish_PS_module(Dict{Symbol, Any}(:figname => pato), cmd, "", true, false, true, G)
+	Ir = gdalread(pato, layout="TCBa")
 	nr_r, nc_r = size(Ir,1), size(Ir,2)
 
 	center_left  = (nr_l ÷ 2 + 1, nc_l ÷ 2 + 1)

src/extras/isoutlier.jl

@@ -97,8 +97,9 @@ function isoutlier(D::GMTdataset{<:Real,2}; critic=3.0, method::Symbol=:median,
 	(threshold === nothing && width <= 0) && error("Must provide the window length of via the `width` option or use the `threshold` option.")
 	if (width > 0)
 		method in [:median, :mean] || throw(ArgumentError("Unknown method: $method. Here use only one of :median or :mean"))
-		(method == :mean) && (method = boxcar)		# :boxcar is the name in GMT for 'mean'
-		Dres = filter1d(D, filter=(type=method, width=width, highpass=true), E=true)
+		_method = (method == :mean) ? :boxcar : method		# :boxcar is the name in GMT for 'mean'
+		#Dres = filter1d(D, filter=(type=_method, width=width, highpass=true), E=true, Vd=1)
+		Dres = gmt("filter1d -E -F" * string(_method)[1] * "$(width)+h", D)
 		isoutlier(view(Dres.data, :, 2), critic=critic, method=method, threshold=threshold)
 	else
 		isoutlier(view(D.data, :, 2), critic=critic, method=method, threshold=threshold)

src/extras/mapsize2region.jl

@@ -0,0 +1,61 @@
+
+# ---------------------------------------------------------------------------------------------------
+"""
+    region_opt, proj_opt = mapsize2region(proj=?, scale=?, clon=?, clat=?, width=?, height=0, bnds="", plot=false)
+
+Compute the region for a map of user specified projection, scale, width and height.
+
+### Kwargs
+- `proj`: Projection. Either a string (e.g. "m" for Mercator, "lambert", "stere", etc.), or a tuple for
+   projections that demand more parameters. Same syntax as `coast`, etc. 
+- `scale`: Scale of the map in the form of "1:xxxx" (e.g. "1:1000000").
+- `clon`: Center longitude of the map in degrees.
+- `clat`: Center latitude of the map in degrees.
+- `width`: Width of the map in centimeters.
+- `height`: Height of the map in centimeters. If not specified, it defaults to `width`.
+- `bnds`: We need a default limits as first approximation to obtain the seeked region. That normally is the global
+   [-180 180 -90 90]. But that might not work for certain cases as for example the Mercator projection Where no poles
+   are allowed). Thoug the Mercator case is dealt in internally, there are other projections that don't allow global
+   limits. In those cases, you will need to specify `bnds` using the same syntax as for the `region` option of all
+   GMT functions.
+- `plot`: If true, we generate an example map representative of the requested region using `coast` function.
+   Use this only to get an idea of the obtained region. No fine tunings are done here. Default is false.
+
+### Returns
+A tuple with the region in the form of "xmin/xmax/ymin/ymax" and the projection string used.
+Use this ooutput as input on all odules that require a `region` and `projection`.
+
+### Credits
+Stolen fom Tim Hume's idea posted in the GMT forum:
+https://forum.generic-mapping-tools.org/t/script-to-create-a-map-with-defined-width-and-height/5909/17?u=joaquim
+
+### Example
+```julia
+opt_R, opt_J = mapsize2region(proj=(name=:tmerc, center=-177), scale="1:10000000", clon=-177, clat=-21, width=15, height=10)
+
+# And now do a simple coastline map showing the region.
+# Note that the same could be achieved by setting the `plot=true` but this shows more explicitly
+# how to use the mapsize2region() output.
+coast(region=opt_R, proj=opt_J, shore=true, show=true)
+```
+"""
+function mapsize2region(; proj="", scale="", clon=NaN, clat=NaN, width=0, height=0, bnds="", plot=false)::Tuple{String, String}
+	(height == 0) && (height = width)		# If height is not specified, use width
+	@assert proj != "" "Projection must be specified"
+	@assert clon != NaN && clat != NaN "Center longitude and latitude must be specified"
+	@assert width > 0 && height > 0 "Width and height must be positive"
+	@assert contains(scale, ':') "Scale must be in the form of '1:xxxx'"
+	(!isa(proj, StrSymb)) && (proj = parse_J(Dict(:J => proj, :scale => scale), "")[1][4:end];	scale="")	# scale is now in J
+	(bnds != "") && (bnds = parse_R(Dict(:R => bnds), "")[4:end])	# If bnds is not empty, parse it
+	(!isa(scale, StrSymb)) && (scale = parse_Scale(Dict(:S => scale), ""))	# If scale is not a StrSymb, parse it
+	(bnds == "" && proj == "m" || startswith(proj, "merc") || startswith(proj, "Merc")) && (bnds="-180/180/-85/85")	# Default bounds for Mercator
+	opt_R, opt_J = mapsize2region(string(proj), scale, Float64(clon), Float64(clat), Float64(width), Float64(height), bnds)
+	(plot != 0) && coast(R=opt_R, J=opt_J, shore=true, show=true, Vd=1)
+	return opt_R, opt_J
+end
+function mapsize2region(proj::String, scale::String, clon::Float64, clat::Float64, width::Float64, height::Float64, bnds)::Tuple{String, String}
+	t::Matrix{Float64} = mapproject([clon clat], J=proj*scale, R=bnds).data
+	t2 = [t[1]-width/2 t[2]-height/2; t[1]+width/2 t[2]+height/2]
+	t3 = mapproject(t2, J=proj*scale, R=bnds, I=true)		# Convert back to lon lat
+	@sprintf("%.8f/%.8f/%.8f/%.8f+r", t3[1,1], t3[1,2], t3[2,1], t3[2,2]), proj*scale
+end

src/gdal_tools.jl

@@ -251,7 +251,7 @@ function gdaldem(indata, method::String, opts::Vector{String}=String[]; dest="/v
 		append!(opts, ["-compute_edges", "-b", band])
 		if ((val = find_in_dict(d, [:scale])[1]) === nothing)
 			if (isa(indata, GMTgrid) && (occursin("longlat", indata.proj4) || occursin("latlong", indata.proj4)) ||
-			                             grdinfo(indata, C="n").data[end] == 1)
+			                             gmt("grdinfo -Cn", indata).data[end] == 1)
 				append!(opts, ["-s", "111120"])
 			end
 		else

src/grdimage.jl

@@ -217,10 +217,6 @@ function common_insert_R!(d::Dict, O::Bool, cmd0, I_G; is3D=false)
 	(opt_R != "") && (CTRL.pocket_R[1] = " -R" * opt_R)
 	return nothing
 end
-function isimgsize(GI)::Bool
-	width, height = getsize(GI)
-	(GI.range[2] - GI.range[1]) == width && (GI.range[4] - GI.range[3]) == height
-end
 
 # ---------------------------------------------------------------------------------------------------
 function common_shade(d::Dict, cmd::String, arg1, arg2, arg3, arg4, prog)

src/histo_funs.jl

@@ -0,0 +1,157 @@
+function find_histo_limits(In, thresholds=nothing, width=20, hst_::Matrix{Float64}=Matrix{Float64}(undef,0,2))
+	# Find the histogram limits of a UInt16 GMTimage that allow to better stretch the histogram
+	# THRESHOLDS is an optional Tuple input containing the left and right histo thresholds, in percentage,
+	# between which the histogram values will be retained. Defaults are (0.1, 0.4) percent. Note, this
+	# assumes the histogram follows some sort of Gaussian distribution. It it's flat, shit occurs.
+	# WIDTH is bin width used to obtain a rough histogram that is used to compute the limits.
+	if (isa(In, Array{UInt16,3}) || isa(In, Array{UInt8,3}))
+		L1 = find_histo_limits(view(In, :, :, 1), thresholds, width)
+		L2 = find_histo_limits(view(In, :, :, 2), thresholds, width)
+		L3 = find_histo_limits(view(In, :, :, 3), thresholds, width)
+		return (L1[1], L1[2], L2[1], L2[2], L3[1], L3[2])
+	end
+	hst = (isempty(hst_)) ? loc_histo(In, "", string(width), "")[1] : hst_
+	if (size(hst, 1) > 10)
+		all(hst[2:5,2] .== 0) && (hst[1,2] = 0)	# Here we always check for high counts in zero bin
+		# Some processed bands leave garbage on the low DNs and that fouls our detecting algo. So check more
+		((hst[1,2] != 0) && hst[1,2] > 100 * mean(hst[2:10,2])) && (hst[1,2] = 0)	# Ad-hoc numbers
+		# Next is for the case of VIIRS bands that have nodata = 65535 but when called from 'truecolor'
+		# it only passes a band view (the array) and we have access to image's header here.
+		(width == 20 && eltype(In) == UInt16 && size(hst,1) == 3277) && (hst[end, 2] = 0)
+	end
+	max_ = maximum(hst, dims=1)[2]
+	(max_ == 0) && error("This histogram had nothing but countings ONLY in first bin. No point to proceed.")
+	thresh_l::Float64 = 0.001;		thresh_r::Float64 = 0.004
+	if (isa(thresholds, Tuple) && length(thresholds) == 2)
+		thresh_l, thresh_r = thresholds[:] ./ 100
+	end
+	thresh_l *= max_
+	thresh_r *= max_
+	kl = 1;		kr = size(hst, 1)
+	while (hst[kl,2] == 0 || hst[kl,2] < thresh_l)  kl += 1  end
+	while (hst[kr,2] == 0 || hst[kr,2] < thresh_r)  kr -= 1  end
+	#return Int(hst[kl,1]), Int(min(hst[kr,1] + width, hst[end,1]))
+	return hst[kl,1], min(hst[kr,1] + width, hst[end,1])
+end
+
+# ---------------------------------------------------------------------------------------------------
+function loc_histo(in, cmd::String="", opt_T::String="", opt_Z::String="")
+	# Very simple function to compute histograms of images (integers)
+	# We put the countings in a Mx2 arrray to trick GMT (pshistogram) to think it's recieving a weighted input.
+	(!isa(in[1], UInt16) && !isa(in[1], UInt8)) && error("Only UInt8 or UInt16 image types allowed here")
+
+	inc::Float64 = (opt_T != "") ? parse(Float64, opt_T) : 1.0
+	(inc <= 0) && error("Bin width must be a number > 0 and no min/max")
+
+	n_bins::Int = (isa(in[1], UInt8)) ? 256 : Int(ceil((maximum(in) + 1) / inc))	# For UInt8 use the full [0 255] range
+	hst = zeros(n_bins, 2)
+	pshst_wall!(in, hst, inc, n_bins)
+
+	cmd = (opt_Z == "") ? cmd * " -Z0" : cmd * opt_Z
+	(!occursin("+w", cmd)) && (cmd *= "+w")			# Pretending to be weighted is crutial for the trick
+
+	return hst, cmd * " -T0/$(n_bins * inc)/$inc"
+end
+
+# ---------------------------------------------------------------------------------------------------
+function pshst_wall!(in, hst, inc, n_bins::Int)
+	# Function barrier for type instability. With the body of this in the calling fun the 'inc' var
+	# introduces a mysterious type instability and execution times multiply by 3.
+	if (inc == 1)
+		@inbounds Threads.@threads for k = 1:numel(in)  hst[in[k] + 1, 2] += 1  end
+    else
+		@inbounds Threads.@threads for k = 1:numel(in)  hst[Int(floor(in[k] / inc) + 1), 2] += 1  end
+	end
+	(isa(in, GItype) && in.nodata == typemax(eltype(in))) && (hst[end] = 0)
+	@inbounds Threads.@threads for k = 1:n_bins  hst[k,1] = inc * (k - 1)  end
+	return nothing
+end
+
+# ---------------------------------------------------------------------------------------------------
+# This version computes the histogram for a UInt8 image band with a bin width of 1
+histogray(img::GMTimage{<:UInt8}; band=1) = histogray(view(img.image, :, :, band))
+function histogray(img::AbstractMatrix{UInt8})
+	edges, counts = 0:255, fill(0, 256)
+	Threads.@threads for v in img
+		@inbounds counts[v+1] += 1
+	end
+	return counts, edges
+end
+
+# ---------------------------------------------------------------------------------------------------
+function hst_floats(arg1, opt_T::String=""; min_max=(0.0, 0.0))
+	# Compute the histogram of a grid or matrix
+	# Made a separate function to let it be called from rescale() and thus avoid calling the main histogram()
+	# that seems to be be too havy (at least according to JET)
+	_min_max::Tuple{Float64,Float64} = (isa(arg1, GMTgrid)) ?
+	                                   (arg1.range[5], arg1.range[6]) : (min_max != (0.0, 0.0) ? min_max : Float64.(extrema_nan(arg1))) 
+	if (opt_T != "")
+		inc = parse(Float64, opt_T) + eps()		# + EPS to avoid the extra last bin at right with 1 count only
+		n_bins = Int(ceil((_min_max[2] - _min_max[1]) / inc))
+	else
+		n_bins = Int(ceil(sqrt(length(arg1))))
+		reg = isa(arg1, GMTgrid) ? (1 - arg1.registration) : 1	# When called from RemoteS arg1 is a view of a layer.
+		inc = (_min_max[2] - _min_max[1]) / (n_bins - reg) + eps()
+	end
+	(!isa(inc, Real) || inc <= 0) && error("Bin width must be a > 0 number and no min/max")
+	hst = zeros(n_bins, 2)
+	have_nans = false
+	if (eltype(arg1) <: AbstractFloat)		# Float arrays can have NaNs
+		have_nans = !(isa(arg1, GMTgrid) && arg1.hasnans == 1)
+		have_nans && (have_nans = any(!isfinite, arg1))
+	end
+
+	_inc = inc + 10eps()					# To avoid cases when index computing fall of by 1
+	if (have_nans)							# If we have NaNs in the grid, we need to take a slower branch
+		@inbounds for k = 1:numel(arg1)
+			!isnan(arg1[k]) && (hst[Int(floor((arg1[k] - _min_max[1]) / _inc) + 1), 2] += 1)
+		end
+	else
+		@inbounds for k = 1:numel(arg1)  hst[Int(floor((arg1[k] - _min_max[1]) / _inc) + 1), 2] += 1  end
+	end
+	@inbounds for k = 1:n_bins  hst[k,1] = _min_max[1] + inc * (k - 1)  end
+	return hst, inc, _min_max
+end
+
+# ---------------------------------------------------------------------------------------------------
+function binmethod(d::Dict, cmd::String, X, is_datetime::Bool)
+	# Compute bin width for a series of binning alghoritms or intervals when X (DateTime) comes in seconds
+	val::String = ((val_ = find_in_dict(d, [:binmethod :BinMethod])[1]) !== nothing) ? lowercase(string(val_)) : ""
+	min_max = (zero(eltype(X)), zero(eltype(X)))
+	(!is_datetime) && (min_max = extrema(X))		# X should already be sorted but don't trust that
+	if (val == "")
+		if (!is_datetime)
+			val = "sqrt"
+		else
+			min_max = extrema(X)		# X should already be sorted but don't trust that
+			rng = (min_max[2] - min_max[1])
+			if     (rng < 150)                 val = "second"
+			elseif (rng / (60) < 150)          val = "minute"
+			elseif (rng / (3600)  < 150)       val = "hour"
+			elseif (rng / (86400) < 150)       val = "day"
+			elseif (rng / (86400 * 7)  < 150)  val = "week"
+			elseif (rng / (86400 * 31) < 150)  val = "month"
+			else                               val = "year"
+			end
+		end
+	end
+
+	n_bins = 0.0;	bin = 0
+	if     (val == "scott")   n_bins = 3.5 .* std(X) .* length(X)^(-1/3)
+	elseif (val == "fd")      n_bins = 2 .* IQR(X) .* length(X)^(-1/3)
+	elseif (val == "sturges") n_bins = ceil.(1 .+ log2.(length(X)))
+	elseif (val == "sqrt")    n_bins = ceil.(sqrt(length(X)))
+	elseif (val == "year")    bin = 86400 * 365.25
+	elseif (val == "month")   bin = 86400 * 31
+	elseif (val == "week")    bin = 86400 * 7
+	elseif (val == "day")     bin = 86400
+	elseif (val == "hour")    bin = 3600
+	elseif (val == "minute")  bin = 60
+	elseif (val == "second")  bin = 1
+	elseif (!is_datetime)     error("Unknown BinMethod $val")
+	end
+	if (bin == 0)
+		bin = (min_max[2] - min_max[1]) / n_bins	# Should be made a "pretty" number?
+	end
+	return @sprintf("%.12g", bin), min_max
+end

src/lepto_funs.jl

@@ -368,11 +368,13 @@ function fillsinks(G::GMTgrid; conn=4, region=nothing, saco=false, insitu=false)
 	d = (I .== I2')
 	D = polygonize(d)
 	if (saco == 1)
-		Dtrk = grdtrack(G, D)
+		Dtrk = gmt("grdtrack -G -n+a", D, G)
 		means = isa(D, Vector) ? median.(Dtrk) : median(Dtrk)
 		SACO[1] = Dict("saco" => Dtrk)			# Save the grdtrack interpolated lines for eventual external use.
 	else
-		means = isa(D, Vector) ? median.(grdtrack(G, D, o=2)) : median(grdtrack(G, D, o=2))	# The mean of each interpolated contour
+		#Dtrk = grdtrack(G, D, o=2)
+		Dtrk = gmt("grdtrack -G -n+a -o2", D, G)
+		means = isa(D, Vector) ? median.(Dtrk) : median(Dtrk)	# The mean of each interpolated contour
 	end
 	_G = (insitu == 1) ? G : deepcopy(G)