Let the getbyattrib() also filter by the largest polygon of each group. (#1789)

* Let gmt_centroid_area select if return centroid,area or just one of them

* Let the getbyattrib() also filter by the largest polygon of each group.

* Add a scatter method to let do choropleths by painted symbols

Author: joa-quim

src/libgmt.jl

@@ -427,31 +427,39 @@ 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}
+function gmt_centroid_area(API::Ptr{Cvoid}, D, geo::Int=0; ca::Int=3)::Matrix{Float64}
+	# geo = 0 for Cartesian, !=0 for geographic
+	# ca = 3 for (x,y,area), 2 for (x,y), 1 for area only
+	mat = Matrix{Float64}(undef, length(D), ca)
 	if (isa(D, Vector))		# 'D' is a vector of GMTdatasets but that type is not yet known here.
-		mat = Matrix{Float64}(undef, length(D), 3)
 		for k in eachindex(D)
 			ref1, ref2, nr = pointercols(D[k])
 			t = gmt_centroid_area(API, ref1, ref2, geo, nr)
-			mat[k, 1], mat[k, 2], mat[k, 3] = t[1], t[2], t[3]
+			if (ca == 3)      mat[k, 1], mat[k, 2], mat[k, 3] = t[1], t[2], t[3]	# To avoid allocs
+			elseif (ca == 2)  mat[k, 1], mat[k, 2] = t[1], t[2]
+			else              mat[k] = t[3]
+			end
 		end
-		mat
 	else
 		ref1, ref2, nr = pointercols(D)
-		gmt_centroid_area(API, ref1, ref2, geo, nr)
+		t = gmt_centroid_area(API, ref1, ref2, geo, nr)[inds]
+		mat = bla(mat, t, ca, 1)
+		if (ca == 3) mat[1], mat[2], mat[3] = t[1], t[2], t[3] elseif (ca == 2) mat[1], mat[2] = t[1], t[2] else mat[1] = t[3] end
 	end
+	return mat
 end
 
-function gmt_centroid_area(API::Ptr{Cvoid}, x, y, geo::Int=0, n::Int=0)::Matrix{Float64}
+function gmt_centroid_area(API::Ptr{Cvoid}, x, y, geo::Int=0, n::Int=0; ca::Int=3)::Matrix{Float64}
 	# geo = 0 for Cartesian, !=0 for geographic
+	# ca = 3 for (x,y,area), 2 for (x,y), 1 for area only
 	@assert isa(x, Vector{Float64}) || isa(x, Ptr{Float64}) || isa(x, Ref{Float64}) "Bad type for x"
 	!isa(x, Vector{Float64}) && (n == 0) && error("Must provide 'n' when x is not a Vector")
 	(n == 0) && (n = length(x))
 	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, n, geo, pos)
 	pos[3] = abs(area)
-	return pos
+	return (ca == 3) ? pos : (ca == 2) ? pos[1:1,1:2] : pos[1:1,3:3]
 end
 
 function pointercols(D::AbstractArray{Float64,2}, cols=(1,2))

src/plot.jl

@@ -362,15 +362,13 @@ Parameters
     + **x**, **cross**
     + **y**, **y_dash**
 	
-    and select their sizes with the **markersize** or **size** keyword [default is 8p].
+    and select their sizes with the `markersize` or `size` keyword [default is 8p].
     The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
     points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT=:c,)*	
 - **W** | **pen** | **markeredgecolor** | **mec** :: [Type => Str]
 
     Set pen attributes for lines or the outline of symbols
 - $(opt_savefig)
-
-[`GMT man page`]($(GMTdoc)plot.html)
 """
 function scatter(f::Function, range_x=nothing; first=true, kw...)
 	rang = gen_coords4funs(range_x, "x"; kw...)
@@ -395,6 +393,19 @@ scatter!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2, true),
 bubblechart  = scatter		# Alias that supposedly only plots circles
 bubblechart! = scatter!
 
+function scatter(D::Vector{<:GMTdataset{Float64,2}}; first=true, kw...)
+	d = KW(kw)
+	labels = String[]
+	if ((s_val = hlp_desnany_str(d, [:labels], false)) !== "")
+		ts = fish_attrib_in_str(s_val)
+		labels = [D[k].attrib[ts] for k = 1:length(D)]
+	end
+	Dc = mat2ds(gmt_centroid_area(G_API[1], D, Int(isgeog(D)), ca=2), geom=wkbPoint, text=labels)
+	(is_in_dict(d, [:marker, :Marker, :shape]) === nothing) && (d[:marker] = "circ")
+	(is_in_dict(d, [:ms :markersize :MarkerSize :size]) === nothing) && (d[:ms] = "12p")
+	_common_plot_xyz("", Dc, "bubble", !first, true, false, d)
+end
+
 # ------------------------------------------------------------------------------------------------------
 scatter3(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, mat2ds(arg1), "scatter3",  true, true; kw...)
 scatter3!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, mat2ds(arg1), "scatter3",  false, true; kw...)

src/spatial_funs.jl

@@ -178,9 +178,14 @@ NOTE: Instead of ``getbyattrib`` one case use instead ``filter`` (...,`index=fal
   means select all elements from ``Antioquia`` and ``Caldas`` that have the attribute `feature_id` = 0.
 
   A second set of attributes can be used to select elements by region, number of polygon vertices and area.
-  For that, name the keyword with a leading underscore, e.g. `_region`, `_nps`, `_area`. Their values are
+  For that, name the keyword with a leading underscore, e.g. `_region`, `_nps`, `_area`, `_unique`. Their values are
   passed respectively a 4 elements tuple and numbers. E.g. `_region=(-8.0, -7.0, 37.0, 38.0)`, `_nps=100`, `_area=10`.
-  Areas are in square km when input is in geographic coordinates, otherwise squre data unites.
+  Areas are in square km when input is in geographic coordinates, otherwise square data units. The `_unique` case is
+  a bit special and is meant to be used when more than one polygon share the same attribute value (e.g. countries with islands).
+  In that case, set the value of `_unique` to the name of the attribute that is shared by the polygons (e.g. `_unique="NAME"`).
+  By default (e.g. `_unique=true`), the attribute name is `Feature_ID` which is the one used by GMT when creating unique
+  IDs for polygons read from OGR formats (.shp, .geojson, etc). The uniqueness is determined by selecting the polygon
+  with the largest area.
 
 - `invert, or reverse, or not`: If `true` return all segments that do NOT match the query condition(s).
 
@@ -208,7 +213,7 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 	invert = (find_in_kwargs(kw, [:invert :not :revert :reverse])[1] !== nothing)
 	if ((_att = find_in_kwargs(kw, [:att :attrib])[1]) !== nothing)		# For backward compat.
 		if !isa(_att, NamedTuple)
-			((val = find_in_kwargs(kw, [:val :value])[1])  === nothing) && error("Must provide the `attribute` VALUE.")
+			((val = find_in_kwargs(kw, [:val :value])[1]) === nothing) && error("Must provide the `attribute` VALUE.")
 		end
 		atts, vals = isa(_att, NamedTuple) ? (string.(keys(_att)), string.(values(_att))) : ((string(_att),), (string(val),))
 		_keys = repeat(" ", length(kw))
@@ -224,6 +229,7 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 		end
 		atts, vals = Vector{String}(undef, count), Vector{String}(undef, count)
 
+		attrib_name = "Feature_ID"			# The default name for the _unique attribute name
 		for k = 1:numel(kw)
 			vv = values(v[k])
 			if (isa(vv, Tuple) && (_keys[k][1] != '_'))
@@ -236,6 +242,11 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 				kk += length(vv) 
 			else
 				atts[kk], vals[kk] = _keys[k], string(vv)
+				if (atts[kk] == "_unique")
+					# If we can't parse the arg as a number then it must be an attribute name. If not, an error will be raised later.
+					(vals[kk] != "true" && (tryparse(Int, vals[kk]) === nothing) && (tryparse(Float64, vals[kk]) === nothing)) && (attrib_name = vals[kk])
+					vals[kk] = "0"	# If _unique the value doesn't matter but must be parseable to float
+				end
 				kk += 1
 			end
 			(k == 1) && (dounion = kk)		# Index that separates the unions from the interceptions
@@ -256,14 +267,30 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 		for k = 1:numel(D)  _tf[k] = size(D[k].data,1) >= nps_  end
 		_tf
 	end
-	function clip_area(D, areas_, area_, _tf)				# Clip by number of points
+	function clip_area(D, areas_, area_, _tf)	# Clip by number of points
 		for k = 1:numel(D)  _tf[k] = areas_[k] >= area_  end
 		_tf
 	end
+	function clip_unique(D, areas_, _tf, name)	# Clip by uniqueness by using the areas to select the largest by group.
+		att_tbl, att_names = make_attrtbl(D, true)
+		((ind_name = findfirst(att_names .== name)) === nothing) && error("Attribute name $name not found in dataset.")
+		_, ind = gunique(att_tbl[:,ind_name])
+		for k = 1:numel(ind)  _tf[ind[k]] = true  end	# Set the unique values to true.
+		k = 1
+		while (k <= numel(ind)-1)
+			((ind[k+=1] - ind[k-1]) == 1) && continue
+			n_start, n_end = ind[k - 1], ind[k] - 1
+			# Here we have a group (from n_start to n_end) of the same type of attribute
+			this_ind = argmax(view(areas_,n_start:n_end,1)) + n_start - 1	# Get the index of the largest area
+			_tf[n_start], _tf[this_ind] = false, true				# Turn off the first and turn on the largest of this group
+		end
+		_tf
+	end
 
 	(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 = gmt_centroid_area(G_API[1], D, Int(isgeog(D)))[:,3])
+	((ind = findfirst(atts .== "_area"))  !== nothing || (ind = findfirst(atts .== "_unique")) !== nothing) && 
+		(area = parse.(Float64, vals[ind]); areas = gmt_centroid_area(G_API[1], D, Int(isgeog(D)), ca=1); isgeog(D) && (areas .*= 1e-6))	# areas in km^2 if geographic coords
 
 	indices::Vector{Int} = Int[]
 	ky = keys(D[1].attrib)
@@ -274,6 +301,7 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 		if (special)
 			if     (atts[n] == "_region") tf = clip_region(D, lims, tf)
 			elseif (atts[n] == "_area")   tf = clip_area(D, areas, area, tf)
+			elseif (atts[n] == "_unique") tf = clip_unique(D, areas, tf, attrib_name)
 			else                          tf = clip_np(D, nps, tf)
 			end
 		else
@@ -301,7 +329,18 @@ function getbyattrib(D::Vector{<:GMTdataset}; indices=false, kw...)::Union{Nothi
 end
 
 # ---------------------------------------------------------------------------------------------------
+"""
+    filter(D::Vector{<:GMTdataset}; kw...)
+
+See `? getbyattrib`
+"""
 Base.:filter(D::Vector{<:GMTdataset}; kw...) = getbyattrib(D; kw...)
+# ---------------------------------------------------------------------------------------------------
+"""
+    findall(D::Vector{<:GMTdataset}; kw...)
+
+See `? getbyattrib`
+"""
 Base.:findall(D::Vector{<:GMTdataset}; kw...) = getbyattrib(D, true; kw...)
 
 # ---------------------------------------------------------------------------------------------------

src/utils_project.jl

@@ -44,7 +44,7 @@ end
 
 # -----------------------------------------------------------------------------------------------
 """
-    GI[,coast] = worldrectangular(GI; proj::String="+proj=vandg", pm=0, latlim=:auto, coast=false)
+    GI[,coast] = worldrectangular(fname::String|GI::GItype; proj::String="+proj=vandg", pm=0, latlim=:auto, coast=false)
 
 Try to create a rectangular map out miscellaneous and not cylindrical projections.
 
@@ -60,7 +60,7 @@ Try to create a rectangular map out miscellaneous and not cylindrical projection
    other than the GSHHG GMT database.
 
 ### Returns
-A grid or an image and optionally the coastlines ... or errors. Not many projections support the procedure
+A grid or an image and optionally the coastlines. 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
 
@@ -157,12 +157,12 @@ end
 
 # -----------------------------------------------------------------------------------------------
 """
-    GI[,coast] = leepacific(fname::String; latlims=nothing, lonlims=nothing, coast=true)
+    GI[,coast] = leepacific(fname::String|GI::GItype; latlims=nothing, lonlims=nothing, coast=true)
 
 Project a geographical grid/image in the Lee Oblated Stereographic projection centered on the Pacific Ocean.
 
-- `GI`: A GMTgrid or GMTimage data type. `GI` can also be a string with a file name of a grid or image.
-   NOTE: This grid/image should have longitudes covering the range 90 to 300 degrees (here lon in [0 360] is better)
+- `GI/fname`: A GMTgrid or GMTimage data type. `GI` can also be a string with a file name of a grid or image.
+   NOTE: This grid/image should (ideally) have longitudes covering the range 76 to 319 degrees (here lon in [0 360] is better)
 - `latlims`: Latitudes used in `region` when reading the grid/image. The default is (-87, 75).
 - `lonlims`: Longitudes used in `region` when reading the grid/image. You cannot deviate mutch from (90,300).
    Actually, while a bug in GDAL is not fixed, you should not change the default values.
@@ -171,13 +171,17 @@ Project a geographical grid/image in the Lee Oblated Stereographic projection ce
    Pass `coast=D`, where `D` is vector of GMTdataset containing coastline polygons with a provenience
    other than the GSHHG GMT database. If `coast=false` the funtion returns only the projected grid/image.
 
+Note: This function uses ``worldrectangular`` with the `+proj=lee_os` projection. If one wants to add
+other elements to the map (points, lines, etc) one must projectem first with the `+proj=lee_os` projection
+(use ``lonlat2xy`` or ``mapproject`` to that purpose).
+
 ### Returns
 A grid or an image and optionally the coastlines.
 
 ### Example:
    G,cl = leepacific("@earth_relief_10m_g");
    grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none)
-   plotgrid!(G_, show=true)
+   plotgrid!(G, show=true)
 """
 function leepacific(fname::String; region=(90.0, 300.0, -90.0, 75.0), latlims=nothing, lonlims=nothing, coast=true)
 	reg = Float64.(collect(region))