From 6395016d3ae2ee4a1e341717755c68d4e84e3fc8 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 4 Oct 2025 01:28:52 +0100 Subject: [PATCH 1/4] Split of mapsize2region() function --- src/extras/mapsize2region.jl | 61 ++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 src/extras/mapsize2region.jl diff --git a/src/extras/mapsize2region.jl b/src/extras/mapsize2region.jl new file mode 100644 index 000000000..76e9d4735 --- /dev/null +++ b/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 From 9cd08cd41efab1b294e4affd7ace454865a60ec7 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 4 Oct 2025 01:31:32 +0100 Subject: [PATCH 2/4] Split of some functions from the pshistogram.jl file --- src/histo_funs.jl | 157 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 src/histo_funs.jl diff --git a/src/histo_funs.jl b/src/histo_funs.jl new file mode 100644 index 000000000..e0ccd1290 --- /dev/null +++ b/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 From 2ab7169ae0b8f909e010e7d4d6ee64c3dc3163cc Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 4 Oct 2025 01:32:40 +0100 Subject: [PATCH 3/4] Move some functions around and replace some calls to module names by GMT.jl(...) --- src/GMT.jl | 3 + src/beziers.jl | 6 +- src/common_options.jl | 6 +- src/extras/anaglyph.jl | 16 ++-- src/extras/isoutlier.jl | 5 +- src/gdal_tools.jl | 2 +- src/grdimage.jl | 4 - src/lepto_funs.jl | 6 +- src/proj_utils.jl | 4 +- src/pshistogram.jl | 159 ----------------------------------- src/psxy.jl | 177 --------------------------------------- src/spatial_funs.jl | 3 +- src/utils.jl | 69 ++------------- test/runtests.jl | 4 + test/test_B-GMTs.jl | 2 +- test/test_avatars.jl | 3 + test/test_common_opts.jl | 2 +- test/test_utils.jl | 2 - 18 files changed, 50 insertions(+), 423 deletions(-) diff --git a/src/GMT.jl b/src/GMT.jl index 43d0845f5..133fffbbe 100644 --- a/src/GMT.jl +++ b/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") diff --git a/src/beziers.jl b/src/beziers.jl index f63f5c7f4..a41e600be 100644 --- a/src/beziers.jl +++ b/src/beziers.jl @@ -62,7 +62,8 @@ 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 @@ -70,7 +71,8 @@ function bezier(D::GMTdataset; t=nothing, np::Int=0, pure=false, firstcurve=true @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 diff --git a/src/common_options.jl b/src/common_options.jl index fd54171b8..8cc00ac42 100644 --- a/src/common_options.jl +++ b/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 diff --git a/src/extras/anaglyph.jl b/src/extras/anaglyph.jl index ec55d79ce..fa02b09a3 100644 --- a/src/extras/anaglyph.jl +++ b/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) diff --git a/src/extras/isoutlier.jl b/src/extras/isoutlier.jl index a808380c7..7a0ca6f9e 100644 --- a/src/extras/isoutlier.jl +++ b/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) diff --git a/src/gdal_tools.jl b/src/gdal_tools.jl index 381bbab70..b7fa12435 100644 --- a/src/gdal_tools.jl +++ b/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 diff --git a/src/grdimage.jl b/src/grdimage.jl index 52b686e3d..bfcae77a7 100644 --- a/src/grdimage.jl +++ b/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) diff --git a/src/lepto_funs.jl b/src/lepto_funs.jl index 87c8242b4..678fc0852 100644 --- a/src/lepto_funs.jl +++ b/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) diff --git a/src/proj_utils.jl b/src/proj_utils.jl index 1e87d7b05..bc396abcf 100644 --- a/src/proj_utils.jl +++ b/src/proj_utils.jl @@ -1017,9 +1017,9 @@ function helper_geoglimits(prj::String, region::Vector{Float64})::Vector{Float64 Gdal.CPLPushErrorHandler(@cfunction(Gdal.CPLQuietErrorHandler, Cvoid, (UInt32, Cint, Cstring))) opts = ["-s_srs", prj, "-t_srs", prj4WGS84, "-overwrite"] ds = Gdal.get_gdaldataset([region[1] region[3]], opts, false)[1] - o1 = Gdal.gdalvectortranslate(ds, opts; dest="/vsimem/tmp", gdataset=true)::GMT.Gdal.IDataset + o1 = Gdal.gdalvectortranslate(ds, opts; dest="/vsimem/tmp", gdataset=true)::Gdal.IDataset ds = Gdal.get_gdaldataset([region[2] region[4]], opts, false)[1] - o2 = Gdal.gdalvectortranslate(ds, opts; dest="/vsimem/tmp", gdataset=true)::GMT.Gdal.IDataset + o2 = Gdal.gdalvectortranslate(ds, opts; dest="/vsimem/tmp", gdataset=true)::Gdal.IDataset if (o1.ptr == C_NULL || o2.ptr == C_NULL) # Diagonals failed, probably a Mollweide or alike projection t = xy2lonlat([region[1] min(0,region[4]); region[2] min(0,region[4]); 0 region; 0 region], s_srs=prj, t_srs=prj4WGS84) diff --git a/src/pshistogram.jl b/src/pshistogram.jl index 9d9b4ad43..2ae6cbfa8 100644 --- a/src/pshistogram.jl +++ b/src/pshistogram.jl @@ -280,165 +280,6 @@ function three_histos(d::Dict, I::GMTimage{UInt8, 3}, cmd, proggy, O, opt_T, opt prep_and_call_finish_PS_module(d, [_cmd * " -Y$(H)c"], "", true, O, true, hst) end -# --------------------------------------------------------------------------------------------------- -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 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 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 - -# --------------------------------------------------------------------------------------------------- -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 - -# --------------------------------------------------------------------------------------------------- -# 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 - # --------------------------------------------------------------------------------------------------- const pshistogram = histogram # Alias const pshistogram! = histogram! # Alias \ No newline at end of file diff --git a/src/psxy.jl b/src/psxy.jl index c6d36d47b..52c804428 100644 --- a/src/psxy.jl +++ b/src/psxy.jl @@ -1350,183 +1350,6 @@ function make_color_column_(d::Dict, cmd::String, len_cmd::Int, N_args::Int, n_p return cmd, arg1, arg2, N_args, true end -# --------------------------------------------------------------------------------------------------- -function get_marker_name(d::Dict, arg1, symbs::Vector{Symbol}, is3D::Bool, del::Bool=true) - marca::String = ""; N = 0 - for symb in symbs - if (haskey(d, symb)) - t = d[symb] - if (isa(t, Tuple)) # e.g. marker=(:r, [2 3]) - msg = ""; cst = false - opt = "" # Probably this defaut value is never used but avoids compiling helper_markers(opt,) with a non def var - o::String = string(t[1]) - if (startswith(o, "E")) opt = "E"; N = 3; cst = true - elseif (startswith(o, "e")) opt = "e"; N = 3 - elseif (o == "J" || startswith(o, "Rot")) opt = "J"; N = 3; cst = true - elseif (o == "j" || startswith(o, "rot")) opt = "j"; N = 3 - elseif (o == "M" || startswith(o, "Mat")) opt = "M"; N = 3 - elseif (o == "m" || startswith(o, "mat")) opt = "m"; N = 3 - elseif (o == "R" || startswith(o, "Rec")) opt = "R"; N = 3 - elseif (o == "r" || startswith(o, "rec")) opt = "r"; N = 2 - elseif (o == "V" || startswith(o, "Vec")) opt = "V"; N = 2 - elseif (o == "v" || startswith(o, "vec")) opt = "v"; N = 2 - elseif (startswith(o, "geovec")) opt = "="; N = 2 - elseif (o == "w" || o == "pie" || o == "web" || o == "wedge") opt = "w"; N = 2 - elseif (o == "W" || o == "Pie" || o == "Web" || o == "Wedge") opt = "W"; N = 2 - end - if (N > 0) marca, arg1, msg = helper_markers(opt, t[2], arg1, N, cst) end - (msg != "") && error(msg) - if (length(t) == 3 && isa(t[3], NamedTuple)) - if (marca == "w" || marca == "W") # Ex (spiderweb): marker=(:pie, [...], (inner=1,)) - marca *= add_opt(t[3], (inner="/", arc="+a", radial="+r", size=("", arg2str, 1), pen=("+p", add_opt_pen)) ) - elseif (marca == "m" || marca == "M" || marca == "=") - marca *= vector_attrib(t[3]) - end - end - elseif (isa(t, NamedTuple)) # e.g. marker=(pie=true, inner=1, ...) - key = keys(t)[1]; opt = "" - if (key == :w || key == :pie || key == :web || key == :wedge) opt = "w" - elseif (key == :W || key == :Pie || key == :Web || key == :Wedge) opt = "W" - elseif (key == :b || key == :bar) opt = "b" - elseif (key == :B || key == :HBar) opt = "B" - elseif (key == :l || key == :letter) opt = "l" - elseif (key == :K || key == :Custom) opt = "K" - elseif (key == :k || key == :custom) opt = "k" - elseif (key == :M || key == :Matang) opt = "M" - elseif (key == :m || key == :matang) opt = "m" - elseif (key == :geovec) opt = "=" - end - if (opt == "w" || opt == "W") - marca = opt * add_opt(t, (size=("", arg2str, 1), inner="/", arc="+a", radial="+r", pen=("+p", add_opt_pen))) - elseif (opt == "b" || opt == "B") - marca = opt * add_opt(t, (size=("", arg2str, 1), base="+b", Base="+B")) - elseif (opt == "l") - marca = opt * add_opt(t, (size=("", arg2str, 1), letter="+t", justify="+j", font=("+f", font))) - elseif (opt == "m" || opt == "M" || opt == "=") - marca = opt * add_opt(t, (size=("", arg2str, 1), arrow=("", vector_attrib))) - elseif (opt == "k" || opt == "K") - marca = opt * add_opt(t, (custom="", size="/")) - end - else - t1::String = string(t) - (t1[1] != 'T') && (t1 = lowercase(t1)) - if (t1 == "-" || t1 == "x-dash") marca = "-" - elseif (t1 == "+" || t1 == "plus") marca = "+" - elseif (t1 == "a" || t1 == "*" || t1 == "star") marca = "a" - elseif (t1 == "k" || t1 == "custom") marca = "k" - elseif (t1 == "x" || t1 == "cross") marca = "x" - elseif (is3D && (t1 == "u" || t1 == "cube")) marca = "u" # Must come before next line - elseif (t1[1] == 'c') marca = "c" - elseif (t1[1] == 'd') marca = "d" # diamond - elseif (t1 == "g" || t1 == "octagon") marca = "g" - elseif (t1[1] == 'h') marca = "h" # hexagon - elseif (t1 == "i" || t1 == "inverted_tri") marca = "i" - elseif (t1[1] == 'l') marca = "l" # letter - elseif (t1 == "n" || t1 == "pentagon") marca = "n" - elseif (t1 == "p" || t1 == "." || t1 == "point") marca = "p" - elseif (t1[1] == 's') marca = "s" # square - elseif (t1[1] == 't' || t1 == "^") marca = "t" # triangle - elseif (t1[1] == 'T') marca = "T" # Triangle - elseif (t1[1] == 'y') marca = "y" # y-dash - elseif (t1[1] == 'f') marca = "f" # for Faults in legend - elseif (t1[1] == 'q') marca = "q" # for Quoted in legend - end - t1 = string(t) # Repeat conversion for the case it was lower-cased above - # Still need to check the simpler forms of these - if (marca == "") marca = helper2_markers(t1, ["e", "ellipse"]) end - if (marca == "") marca = helper2_markers(t1, ["E", "Ellipse"]) end - if (marca == "") marca = helper2_markers(t1, ["j", "rotrect"]) end - if (marca == "") marca = helper2_markers(t1, ["J", "RotRect"]) end - if (marca == "") marca = helper2_markers(t1, ["m", "matangle"]) end - if (marca == "") marca = helper2_markers(t1, ["M", "Matangle"]) end - if (marca == "") marca = helper2_markers(t1, ["r", "rectangle"]) end - if (marca == "") marca = helper2_markers(t1, ["R", "RRectangle"]) end - if (marca == "") marca = helper2_markers(t1, ["v", "vector"]) end - if (marca == "") marca = helper2_markers(t1, ["V", "Vector"]) end - if (marca == "") marca = helper2_markers(t1, ["w", "pie", "web"]) end - if (marca == "") marca = helper2_markers(t1, ["W", "Pie", "Web"]) end - end - (del) && delete!(d, symb) - break - end - end - return marca, arg1, N > 0 -end - -function helper_markers(opt::String, ext, arg1::GMTdataset, N::Int, cst::Bool) - # Helper function to deal with the cases where one sends marker's extra columns via command - # Example that will land and be processed here: marker=(:Ellipse, [30 10 15]) - # N is the number of extra columns - marca::String = ""; msg = "" - if (size(ext,2) == N && arg1 !== nothing) # Here ARG1 is supposed to be a matrix that will be extended. - S = Symbol(opt) - t = arg1.data # Because we need to passa matrix to this method of add_opt() - marca, t = add_opt(add_opt, (Dict(S => (par=ext,)), opt, "", [S]), (par="|",), true, t) - arg1.data = t; add2ds!(arg1) - elseif (cst && length(ext) == 1) - marca = opt * "-" * string(ext)::String - else - msg = string("Wrong number of extra columns for marker (", opt, "). Got ", size(ext,2), " but expected ", N) - end - return marca, arg1, msg -end - -function helper2_markers(opt::String, alias::Vector{String})::String - marca = "" - if (opt == alias[1]) # User used only the one letter syntax - marca = alias[1] - else - for k = 2:length(alias) # Loop because of cases like ["w" "pie" "web"] - o2 = alias[k][1:min(2,length(alias[k]))] # check the first 2 chars and Ro, Rotrect or RotRec are all good - #if (startswith(opt, o2)) marca = alias[1]; break end # Good when, for example, marker=:Pie - if (startswith(opt, o2)) # Good when, for example, marker=:Pie - marca = alias[1]; - (opt[end] == '-') && (marca *= '-') - break - end - end - end - - # If we still have found nothing, assume that OPT is a full GMT opt string (e.g. W/5+a30+r45+p2,red) - (marca == "" && opt[1] == alias[1][1]) && (marca = opt) - return marca -end - -# --------------------------------------------------------------------------------------------------- -function seek_custom_symb(marca::AbstractString, with_k::Bool=false)::String - # If 'marca' is a custom symbol, seek it first in GMT.jl share/custom dir. - # Return the full name of the marker plus extension - # The WITH_K arg is to allow calling this fun with a sym name already prefaced with 'k', or not - (with_k && marca[1] != 'k') && return marca # Not a custom symbol, return what we got. - - function find_this_file(pato, symbname) - for (root, dirs, files) in walkdir(pato) - ind = findfirst(startswith.(files, symbname)) - if (ind !== nothing) return joinpath(root, files[ind]) end - end - return "" - end - - s = split(marca, '/') - ind_s = with_k ? 2 : 1 - symbname = s[1][ind_s:end] - cus_path = joinpath(dirname(pathof(GMT))[1:end-4], "share", "custom") - - fullname = find_this_file(cus_path, symbname) - if (fullname == "") - cus_path2 = joinpath(GMTuserdir[1], "cache_csymb") - cus_path2 = replace(cus_path2, "/" => "\\") # Otherwise it will produce currupted PS - fullname = find_this_file(cus_path2, symbname) - end - - (fullname == "") && return marca # Assume it's a custom symbol from the official GMT collection. - - _siz = split(marca, '/')[2] # The custom symbol size - _marca = (with_k ? "k" : "") * fullname * "/" * _siz - (GMTver <= v"6.4" && (length(_marca) - length(_siz) -2) > 62) && warn("Due to a GMT <= 6.4 limitation the length of full (name+path) custom symbol name cannot be longer than 62 bytes.") - return _marca -end - # --------------------------------------------------------------------------------------------------- function check_caller(d::Dict, cmd::String, opt_S::String, opt_W::String, caller::String, g_bar_fill::Vector{String}, O::Bool)::String # Set sensible defaults for the sub-modules "scatter" & "bar" diff --git a/src/spatial_funs.jl b/src/spatial_funs.jl index bf8e707ff..849b6da41 100644 --- a/src/spatial_funs.jl +++ b/src/spatial_funs.jl @@ -566,7 +566,8 @@ function randinpolygon(Din::GDtype; density=0.1, np::Int=0) x = rand(DT, _np) * _dx .+ x0 y = rand(DT, _np) * _dy .+ y0 end - gmtselect([x y], polygon=_D) + #gmtselect([x y], polygon=_D, Vd=1) + gmt("gmtselect -F", [x y], _D) end isgeo = isgeog(Din) diff --git a/src/utils.jl b/src/utils.jl index e10d48ca6..fee341484 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1203,6 +1203,13 @@ function getface(FV::GMTfv, face=1, n=1; view=false) (view == 1) ? FV.verts[FV.faces_view[face][n,:],:] : FV.verts[FV.faces[face][n,:],:] end +# ------------------------------------------------------------------------------------------------------ +function isimgsize(GI)::Bool + # To find out if coordinates are in fact image sizes, which is used to NOT plot axes when plotting images. + width, height = getsize(GI) + (GI.range[2] - GI.range[1]) == width && (GI.range[4] - GI.range[3]) == height +end + # ------------------------------------------------------------------------------------------------------ """ isclockwise(poly::Matrix{<:AbstractFloat}, view=(0.0,0.0,1.0)) -> Bool @@ -1458,68 +1465,6 @@ function bissextile(year::Integer)::Bool (year % 4 == 0) && ((year % 100 != 0) || (year % 400 == 0)) end - -# --------------------------------------------------------------------------------------------------- -""" - 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 - # --------------------------------------------------------------------------------------------------- include("makeDCWs.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 64e1d3080..6b9655b1b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,6 +60,10 @@ using InteractiveUtils catch err @warn("Failed the WMS test. Error was:\n $err") end + + println(" MAPSIZE2REGION") + # This fck one is extremly irritating. It errors with different errors depending on where is executed. Even with a resetGMT(). + @test mapsize2region(proj=(name=:tmerc, center=-177), scale="1:10000000", clon=-177, clat=-21, width=15, height=10)[2] == "t-177/1:10000000" include("test_isoutlier.jl") include("test_okadas.jl") diff --git a/test/test_B-GMTs.jl b/test/test_B-GMTs.jl index 14189d355..b47b2249b 100644 --- a/test/test_B-GMTs.jl +++ b/test/test_B-GMTs.jl @@ -17,7 +17,7 @@ D = blockmode(region=[0 2 0 2], inc=1, reg=true, d); println(" CONTOURF") - G = GMT.peaks(); + G = peaks(); gmtwrite("lixo.grd", G) C = makecpt(T=(-7,9,2)); contourf("lixo.grd", Vd=dbg2) diff --git a/test/test_avatars.jl b/test/test_avatars.jl index 532e962fe..408a3335f 100644 --- a/test/test_avatars.jl +++ b/test/test_avatars.jl @@ -296,6 +296,7 @@ bar(rand(20),hbar="0.5c+b9", Vd=dbg2) bar(rand(10), xaxis=(custom=(pos=1:5,type="A",angle=45),), Vd=dbg2) bar(rand(10), axis=(custom=(pos=1:5,label=[:a :b :c :d :e]),), Vd=dbg2) + @info "BAR1" bar((1,2,3), Vd=dbg2) bar((1,2,3), (1,2,3), Vd=dbg2) bar!((1,2,3), Vd=dbg2) @@ -305,6 +306,7 @@ bar!("", [3 31], C=:lightblue, frame=:noannot, Vd=dbg2) men_means, men_std = (20, 35, 30, 35, 27), (2, 3, 4, 1, 2); x = collect(1:length(men_means)); + @info "BAR2" bar(x.-0.35/2, collect(men_means), width=0.35, color=:lightblue, limits=(0.5,5.5,0,40), frame=:none, error_bars=(y=men_std,), Vd=dbg2) bar([0. 1 2 3; 1 2 3 4], error_bars=(y=[0.1 0.2 0.33; 0.2 0.3 0.4],), fillalpha=[0.3 0.5 0.7], Vd=dbg2) bar([0. 1 2 3; 1 2 3 4], fill=(1,2,3), Vd=dbg2) @@ -314,6 +316,7 @@ bar(1:3,[-5 -15 20; 17 10 21; 10 5 15], stacked=1, Vd=dbg2) bar([0. 1 2 3; 1 2 3 4], fill=("red", "green", "blue"), Vd=dbg2) bar(rand(15), color=:rainbow, Y=3) + @info "BAR3" T = mat2ds([1.0 0.446143; 2.0 0.581746; 3.0 0.268978], text=[" "; " "; " "]); bar(T, color=:rainbow, figsize=(14,8), title="Colored bars", Vd=dbg2) T = mat2ds([1.0 0.446143 0; 2.0 0.581746 0; 3.0 0.268978 0], text=[" "; " "; " "]); diff --git a/test/test_common_opts.jl b/test/test_common_opts.jl index 9d1e0c974..c1f90a144 100644 --- a/test/test_common_opts.jl +++ b/test/test_common_opts.jl @@ -341,7 +341,7 @@ @test o == " -TG -W+k+tLoLo" o = GMT.prepare2geotif(Dict{Symbol,Any}(:kml => (title=:Lolo, layer=:bla, fade=(1,2), URL="http")), ["pscoast -Rd -JX12 -Baf -W0.5p -Da"], "", false)[2]; @test o == " -TG -W+k+tLolo+nbla+f1/2+uhttp" - coast(region=:global, kml=:trans, proj=:merc,Vd=dbg2) + #coast(region=:global, kml=:trans, proj=:merc,Vd=dbg2) @test (GMT.check_flipaxes(Dict{Symbol,Any}(:flipaxes => (y=1,)), "?") == "-?") @test (GMT.check_flipaxes(Dict{Symbol,Any}(:flipaxes => ("x", :y)), "?/?") == "-?/-?") diff --git a/test/test_utils.jl b/test/test_utils.jl index 9d40111aa..9e2e82d5c 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -42,6 +42,4 @@ @test !bissextile(100) @test bissextile(-4) - - @test mapsize2region(proj=(name=:tmerc, center=-177), scale="1:10000000", clon=-177, clat=-21, width=15, height=10)[2] == "t-177/1:10000000" end From b311e5f9e6f0dcf0fcdbdc18e8fea1452fb9b8b1 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sat, 4 Oct 2025 01:39:40 +0100 Subject: [PATCH 4/4] A split from psxy.jl --- src/marker_name.jl | 176 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 src/marker_name.jl diff --git a/src/marker_name.jl b/src/marker_name.jl new file mode 100644 index 000000000..e91783c06 --- /dev/null +++ b/src/marker_name.jl @@ -0,0 +1,176 @@ +# Put this code in a separate file to permit access from future GMT_base module +function get_marker_name(d::Dict, arg1, symbs::Vector{Symbol}, is3D::Bool, del::Bool=true) + marca::String = ""; N = 0 + for symb in symbs + if (haskey(d, symb)) + t = d[symb] + if (isa(t, Tuple)) # e.g. marker=(:r, [2 3]) + msg = ""; cst = false + opt = "" # Probably this default value is never used but avoids compiling helper_markers(opt,) with a non def var + o::String = string(t[1]) + if (startswith(o, "E")) opt = "E"; N = 3; cst = true + elseif (startswith(o, "e")) opt = "e"; N = 3 + elseif (o == "J" || startswith(o, "Rot")) opt = "J"; N = 3; cst = true + elseif (o == "j" || startswith(o, "rot")) opt = "j"; N = 3 + elseif (o == "M" || startswith(o, "Mat")) opt = "M"; N = 3 + elseif (o == "m" || startswith(o, "mat")) opt = "m"; N = 3 + elseif (o == "R" || startswith(o, "Rec")) opt = "R"; N = 3 + elseif (o == "r" || startswith(o, "rec")) opt = "r"; N = 2 + elseif (o == "V" || startswith(o, "Vec")) opt = "V"; N = 2 + elseif (o == "v" || startswith(o, "vec")) opt = "v"; N = 2 + elseif (startswith(o, "geovec")) opt = "="; N = 2 + elseif (o == "w" || o == "pie" || o == "web" || o == "wedge") opt = "w"; N = 2 + elseif (o == "W" || o == "Pie" || o == "Web" || o == "Wedge") opt = "W"; N = 2 + end + if (N > 0) marca, arg1, msg = helper_markers(opt, t[2], arg1, N, cst) end + (msg != "") && error(msg) + if (length(t) == 3 && isa(t[3], NamedTuple)) + if (marca == "w" || marca == "W") # Ex (spiderweb): marker=(:pie, [...], (inner=1,)) + marca *= add_opt(t[3], (inner="/", arc="+a", radial="+r", size=("", arg2str, 1), pen=("+p", add_opt_pen)) ) + elseif (marca == "m" || marca == "M" || marca == "=") + marca *= vector_attrib(t[3]) + end + end + elseif (isa(t, NamedTuple)) # e.g. marker=(pie=true, inner=1, ...) + key = keys(t)[1]; opt = "" + if (key == :w || key == :pie || key == :web || key == :wedge) opt = "w" + elseif (key == :W || key == :Pie || key == :Web || key == :Wedge) opt = "W" + elseif (key == :b || key == :bar) opt = "b" + elseif (key == :B || key == :HBar) opt = "B" + elseif (key == :l || key == :letter) opt = "l" + elseif (key == :K || key == :Custom) opt = "K" + elseif (key == :k || key == :custom) opt = "k" + elseif (key == :M || key == :Matang) opt = "M" + elseif (key == :m || key == :matang) opt = "m" + elseif (key == :geovec) opt = "=" + end + if (opt == "w" || opt == "W") + marca = opt * add_opt(t, (size=("", arg2str, 1), inner="/", arc="+a", radial="+r", pen=("+p", add_opt_pen))) + elseif (opt == "b" || opt == "B") + marca = opt * add_opt(t, (size=("", arg2str, 1), base="+b", Base="+B")) + elseif (opt == "l") + marca = opt * add_opt(t, (size=("", arg2str, 1), letter="+t", justify="+j", font=("+f", font))) + elseif (opt == "m" || opt == "M" || opt == "=") + marca = opt * add_opt(t, (size=("", arg2str, 1), arrow=("", vector_attrib))) + elseif (opt == "k" || opt == "K") + marca = opt * add_opt(t, (custom="", size="/")) + end + else + t1::String = string(t) + (t1[1] != 'T') && (t1 = lowercase(t1)) + if (t1 == "-" || t1 == "x-dash") marca = "-" + elseif (t1 == "+" || t1 == "plus") marca = "+" + elseif (t1 == "a" || t1 == "*" || t1 == "star") marca = "a" + elseif (t1 == "k" || t1 == "custom") marca = "k" + elseif (t1 == "x" || t1 == "cross") marca = "x" + elseif (is3D && (t1 == "u" || t1 == "cube")) marca = "u" # Must come before next line + elseif (t1[1] == 'c') marca = "c" + elseif (t1[1] == 'd') marca = "d" # diamond + elseif (t1 == "g" || t1 == "octagon") marca = "g" + elseif (t1[1] == 'h') marca = "h" # hexagon + elseif (t1 == "i" || t1 == "inverted_tri") marca = "i" + elseif (t1[1] == 'l') marca = "l" # letter + elseif (t1 == "n" || t1 == "pentagon") marca = "n" + elseif (t1 == "p" || t1 == "." || t1 == "point") marca = "p" + elseif (t1[1] == 's') marca = "s" # square + elseif (t1[1] == 't' || t1 == "^") marca = "t" # triangle + elseif (t1[1] == 'T') marca = "T" # Triangle + elseif (t1[1] == 'y') marca = "y" # y-dash + elseif (t1[1] == 'f') marca = "f" # for Faults in legend + elseif (t1[1] == 'q') marca = "q" # for Quoted in legend + end + t1 = string(t) # Repeat conversion for the case it was lower-cased above + # Still need to check the simpler forms of these + if (marca == "") marca = helper2_markers(t1, ["e", "ellipse"]) end + if (marca == "") marca = helper2_markers(t1, ["E", "Ellipse"]) end + if (marca == "") marca = helper2_markers(t1, ["j", "rotrect"]) end + if (marca == "") marca = helper2_markers(t1, ["J", "RotRect"]) end + if (marca == "") marca = helper2_markers(t1, ["m", "matangle"]) end + if (marca == "") marca = helper2_markers(t1, ["M", "Matangle"]) end + if (marca == "") marca = helper2_markers(t1, ["r", "rectangle"]) end + if (marca == "") marca = helper2_markers(t1, ["R", "RRectangle"]) end + if (marca == "") marca = helper2_markers(t1, ["v", "vector"]) end + if (marca == "") marca = helper2_markers(t1, ["V", "Vector"]) end + if (marca == "") marca = helper2_markers(t1, ["w", "pie", "web"]) end + if (marca == "") marca = helper2_markers(t1, ["W", "Pie", "Web"]) end + end + (del) && delete!(d, symb) + break + end + end + return marca, arg1, N > 0 +end + +function helper_markers(opt::String, ext, arg1::GMTdataset, N::Int, cst::Bool) + # Helper function to deal with the cases where one sends marker's extra columns via command + # Example that will land and be processed here: marker=(:Ellipse, [30 10 15]) + # N is the number of extra columns + marca::String = ""; msg = "" + if (size(ext,2) == N && arg1 !== nothing) # Here ARG1 is supposed to be a matrix that will be extended. + S = Symbol(opt) + t = arg1.data # Because we need to passa matrix to this method of add_opt() + marca, t = add_opt(add_opt, (Dict(S => (par=ext,)), opt, "", [S]), (par="|",), true, t) + arg1.data = t; add2ds!(arg1) + elseif (cst && length(ext) == 1) + marca = opt * "-" * string(ext)::String + else + msg = string("Wrong number of extra columns for marker (", opt, "). Got ", size(ext,2), " but expected ", N) + end + return marca, arg1, msg +end + +function helper2_markers(opt::String, alias::Vector{String})::String + marca = "" + if (opt == alias[1]) # User used only the one letter syntax + marca = alias[1] + else + for k = 2:length(alias) # Loop because of cases like ["w" "pie" "web"] + o2 = alias[k][1:min(2,length(alias[k]))] # check the first 2 chars and Ro, Rotrect or RotRec are all good + #if (startswith(opt, o2)) marca = alias[1]; break end # Good when, for example, marker=:Pie + if (startswith(opt, o2)) # Good when, for example, marker=:Pie + marca = alias[1]; + (opt[end] == '-') && (marca *= '-') + break + end + end + end + + # If we still have found nothing, assume that OPT is a full GMT opt string (e.g. W/5+a30+r45+p2,red) + (marca == "" && opt[1] == alias[1][1]) && (marca = opt) + return marca +end + +# --------------------------------------------------------------------------------------------------- +function seek_custom_symb(marca::AbstractString, with_k::Bool=false)::String + # If 'marca' is a custom symbol, seek it first in GMT.jl share/custom dir. + # Return the full name of the marker plus extension + # The WITH_K arg is to allow calling this fun with a sym name already prefaced with 'k', or not + (with_k && marca[1] != 'k') && return marca # Not a custom symbol, return what we got. + + function find_this_file(pato, symbname) + for (root, dirs, files) in walkdir(pato) + ind = findfirst(startswith.(files, symbname)) + if (ind !== nothing) return joinpath(root, files[ind]) end + end + return "" + end + + s = split(marca, '/') + ind_s = with_k ? 2 : 1 + symbname = s[1][ind_s:end] + cus_path = joinpath(dirname(pathof(GMT))[1:end-4], "share", "custom") + + fullname = find_this_file(cus_path, symbname) + if (fullname == "") + cus_path2 = joinpath(GMTuserdir[1], "cache_csymb") + cus_path2 = replace(cus_path2, "/" => "\\") # Otherwise it will produce currupted PS + fullname = find_this_file(cus_path2, symbname) + end + + (fullname == "") && return marca # Assume it's a custom symbol from the official GMT collection. + + _siz = split(marca, '/')[2] # The custom symbol size + _marca = (with_k ? "k" : "") * fullname * "/" * _siz + (GMTver <= v"6.4" && (length(_marca) - length(_siz) -2) > 62) && warn("Due to a GMT <= 6.4 limitation the length of full (name+path) custom symbol name cannot be longer than 62 bytes.") + return _marca +end