Add a gmt_centroid_area() fun to directly ccall the GMT C lib function. (#1785)

Author: joa-quim

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))...,)))

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

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)

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)
 """

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