diff --git a/src/libgmt.jl b/src/libgmt.jl index 36dde6464..29c54276a 100644 --- a/src/libgmt.jl +++ b/src/libgmt.jl @@ -427,6 +427,28 @@ function gmt_free_mem(API::Ptr{Cvoid}, mem) ccall((:gmt_free_func, libgmt), Cvoid, (Cstring, Ptr{Cvoid}, Bool, Cstring), GMT_, mem, true, "Julia") end +function gmt_centroid_area(API::Ptr{Cvoid}, D, geo::Int=0)::Matrix{Float64} + if (isa(D, Vector)) + mat = Matrix{Float64}(undef, length(D), 3) + for k in eachindex(D) + t = gmt_centroid_area(API::Ptr{Cvoid}, D[k].data[:,1], D[k].data[:,2], geo) + mat[k, 1], mat[k, 2], mat[k, 3] = t[1], t[2], t[3] + end + mat + else + gmt_centroid_area(API::Ptr{Cvoid}, D.data[:,1], D.data[:,2], geo) + end +end + +function gmt_centroid_area(API::Ptr{Cvoid}, x::Vector{Float64}, y::Vector{Float64}, geo::Int=0)::Matrix{Float64} + # geo = 0 for Cartesian, !=0 for geographic + pos = [0.0 0.0 0.0] + area = ccall((:gmt_centroid_area, libgmt), Cdouble, (Cstring, Ptr{Cdouble}, Ptr{Cdouble}, UInt64, Int32, Ptr{Cdouble}), + GMT_Get_Ctrl(API), x, y, length(x), geo, pos) + pos[3] = abs(area) + return pos +end + #= function get_common_R(API::Ptr{Cvoid}) R = COMMON_R((false,false,false,false), false, 0, 0, 0, (0., 0., 0., 0., 0., 0.), (0., 0., 0., 0.), (0., 0.), map(UInt8, (string(repeat(" ",256))...,))) diff --git a/src/psxy.jl b/src/psxy.jl index 3826ac59b..a91f97e6d 100644 --- a/src/psxy.jl +++ b/src/psxy.jl @@ -418,18 +418,19 @@ function plt_txt_attrib!(D::Vector{<:GMTdataset{T,N}}, d::Dict{Symbol, Any}, _cm ((_ind = findfirst('=', s_val)) === nothing) && (_ind = findfirst(':', s_val)) ind::Int = _ind # Because it fck insists _ind is a Any ts::String = s_val[ind+1:end] - ct::GMTdataset = gmtspatial(D, centroid=true) # Texts will be plotted at the polygons centroids + t = vec(make_attrtbl(D, att=ts)[1]) if ((fnt = add_opt(d, "", "", [:font], (angle="+a", font=("+f", font)); del=false)) != "") (fnt[1] != '+') && (fnt = "+f" * fnt) delete!(d, :font) - ct.text = vec(make_attrtbl(D, att=ts)[1]) else nc::Int = round(Int, sqrt(length(D))) # A crude guess of the number of columns fnt = (nc < 5) ? "7p" : (nc < 9 ? "5p" : "4p") # A simple heuristic outline = fnt * ",black=~0.75p,white " # Apply the outline trick fnt = "+f" - ct.text = outline .* vec(make_attrtbl(D, att=ts)[1]) + t = outline .* t end + ct::GMTdataset{Float64,2} = mat2ds(gmt_centroid_area(G_API[1], D, Int(isgeog(D)))) + ct.text = t # Texts will be plotted at the polygons centroids (CTRL.pocket_call[1] === nothing) ? (CTRL.pocket_call[1] = ct) : (CTRL.pocket_call[2] = ct) append!(_cmd, ["pstext -R -J -F" * fnt * "+jMC"]) return nothing diff --git a/src/spatial_funs.jl b/src/spatial_funs.jl index 41df4028d..dd9065924 100644 --- a/src/spatial_funs.jl +++ b/src/spatial_funs.jl @@ -263,7 +263,7 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int} (ind = findfirst(atts .== "_region")) !== nothing && (lims = parse.(Float64, split(vals[ind][2:end-1], ", "))) (ind = findfirst(atts .== "_nps")) !== nothing && (nps = parse.(Float64, vals[ind])) - (ind = findfirst(atts .== "_area")) !== nothing && (area = parse.(Float64, vals[ind]); areas::Vector{Float64} = gmtspatial(D, area=true).data[:,3]) + (ind = findfirst(atts .== "_area")) !== nothing && (area = parse.(Float64, vals[ind]); areas = gmt_centroid_area(G_API[1], D, Int(isgeog(D)))[:,3]) indices::Vector{Int} = Int[] ky = keys(D[1].attrib) diff --git a/src/utils_project.jl b/src/utils_project.jl index 8194ae94f..75c292d96 100644 --- a/src/utils_project.jl +++ b/src/utils_project.jl @@ -171,12 +171,11 @@ Project a geographical grid/image in the Lee Oblated Stereographic projection ce ### Returns -A grid or an image and optionally the coastlines ... or errors. Not many projections support the procedure -implemented in this function. -The working or not is controlled by PROJ's `+over` option https://proj.org/usage/projections.html#longitude-wrapping +A grid or an image and optionally the coastlines. ### Example: - G,cl = leepacific("@earth_relief_10m_g") + G,cl = leepacific("@earth_relief_10m_g"); + grid = worldrectgrid(G); grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none) plotgrid!(G, grid, show=true) """ diff --git a/test/test_new_projs.jl b/test/test_new_projs.jl index 381f78940..7b344cc71 100644 --- a/test/test_new_projs.jl +++ b/test/test_new_projs.jl @@ -8,6 +8,7 @@ G = worldrectangular("@earth_relief_10m_g", pm=80); G = worldrectangular("@earth_relief_10m_g", pm=-100); G,cl = worldrectangular("@earth_relief_10m_g", latlim=90, coast=true); + G,cl = leepacific("@earth_relief_10m_g"); grid = worldrectgrid(G, annot_x=[-180,-150,0,150,180]) plot([0 0]) # Just to have a PS for the next have where to append