Add a new filter option to filter Datasets by BBox aspect ratio.

src/GMT.jl

@@ -196,6 +196,7 @@ export
 	okada,
 	haralick,
 	mapsize2region,
+	ind2bool,
 
 	cube, cylinder, circlepts, dodecahedron, ellipse3D, eulermat, flatfv, icosahedron, loft, sphere, spinmat,
 	octahedron, tetrahedron, torus, replicant, revolve, rotate, rotate!, translate, translate!,

src/imgmorph/cc2bw.jl

@@ -2,20 +2,38 @@
     BW = cc2bw(cc::GMTConComp; obj2keep::Union{Int, Vector{Int}}=0)
 
 Convert connected components to binary image
+
+- `cc`: Connected components created with `bwconncomp`
+- `obj2keep`: (optional) Integer, vector of integers or a vector of booleans specifying which
+   connected components to keep in the output binary image. By default, all components are kept.
+
+Returns a binary `GMTimage` where pixels belonging to the specified connected components are set to `true`.
 """
-function cc2bw(cc::GMTConComp; obj2keep::Union{Int, Vector{Int}}=0)
+function cc2bw(cc::GMTConComp; obj2keep=Int[])
+	isscalar(obj2keep) && return _cc2bw(cc, [Int(obj2keep)])
+	if (isa(obj2keep, Vector{Int}))
+		return _cc2bw(cc, obj2keep)
+	elseif (isa(obj2keep, Vector{Bool}) || isa(obj2keep, BitVector))
+		ind = findall(obj2keep)
+		isempty(ind) && error("No components selected to keep in the 'obj2keep' parameter")
+		return _cc2bw(cc, ind)
+	end
+	error("Invalid type for 'obj2keep' parameter. Must be an integer, vector of integers, or vector of booleans.")
+end
+function _cc2bw(cc::GMTConComp, obj2keep::Vector{Int})
 	# Create a binary image of the same size as the original
+	!isempty(obj2keep) && any((k -> k < 1 || k > cc.num_objects), obj2keep) &&
+	    error("Values in 'obj2keep' are out of range [1, $(cc.num_objects)]")
 	I = mat2img(zeros(Bool, cc.image_size), x=cc.x, y=cc.y, inc=cc.inc, layout=cc.layout,
 	            is_transposed=(cc.layout[2] == 'R'), proj4=cc.proj4, wkt=cc.wkt, epsg=cc.epsg)
 
 	# Mark pixels belonging to any connected component as true
-	plist = (obj2keep == 0) ? (1:cc.num_objects) : isvector(obj2keep) ? obj2keep : (obj2keep:obj2keep)
-	for k = 1:numel(plist)
-		for rc in cc.pixel_list[plist[k]]
+	plist = isempty(obj2keep) ? (1:cc.num_objects) : obj2keep
+	@inbounds for k = 1:numel(plist)
+		@inbounds for rc in cc.pixel_list[plist[k]]
 			I[rc] = true
 		end
 	end
 	I.range[6] = 1.0
 	return I
 end
-

src/lepto_funs.jl

@@ -340,17 +340,6 @@ which is not that much.
 - `conn::Int`: Connectivity for sink filling (4 or 8). Default is 8.
 - `region`: Limit the action to a region of the grid specified by `region`. See for example the ``coast``
   manual for and extended doc on this keword, but note that here only `region` is accepted and not `R`, etc...
-- `saco::Bool`: Save the lines (GMTdataset ~contours) used to fill the sinks in a global variable called
-  SACO. This is intended to avoid return them all the time when function ends. This global variable
-  is a ``[Dict{String,Union{AbstractArray, Vector{AbstractArray}}}()]``, so to access its contents you must use:
-
-  ``D = get(SACO[1], "saco", nothing)``, where ``D`` is now a GMTdataset or a vector of them.
-
-  NOTE: It is the user's responsibility to empty this global variable when it is no longer needed.
-
-  You do that with: ``delete!(SACO[1], "saco")``
-- `insitu::Bool`: If `true`, modify the grid in place. Default is `false`.
-  Alternatively, use the conventional form ``fillsinks!(G; conn=4)``.
 
 ### Returns
 - A new `GMTgrid` with sinks filled, unless `insitu` is `true`, in which case the input grid is modified and returned.
@@ -361,16 +350,8 @@ G = peaks();
 G2 = fillsinks(G);
 viz(G2, shade=true)
 ```
-
-Now save the filling contours and make a plot that overlayes them
-```julia
-G2 = fillsinks(G);
-G2 = fillsinks(G, saco=true);
-grdimage(G2)
-plot!(get(SACO[1], "saco", nothing), color=:white, show=true)
-```
 """
-function fillsinks(GI::GItype; conn=8, region=nothing, saco=false, insitu=false)
+function fillsinks(GI::GItype; conn=8, region=nothing)
 
 	infa = eltype(GI) <: AbstractFloat ? -Inf : typemin(eltype(GI))
 	for k = 1:numel(GI)  GI[k] *= -1  end
@@ -1213,18 +1194,18 @@ A `GMTConComp` structure with the following fields:
 - `connectivity`: Connectivity value used to identify the connected components in I (4 or 8).
 - `num_objects`: Number of connected components found.
 - `image_size`: Size of the image.
-- `range`: Range of the image.
-- `inc`: Increment of the image.
+- `range`: Range of the image coordinates.
+- `inc`: Image's increment (!= 1 whem image is referenced).
 - `registration`: Registration of the image.
 - `x`: X coordinates of the image.
 - `y`: Y coordinates of the image.
-- `layout`: Layout of the image used to find the components.
-- `proj4`: Projection definition of the image used to find the components.
-- `wkt`: Well-known text definition of the image used to find the components.
-- `epsg`: EPSG code of the image used to find the components.
+- `layout`: Memory layout of the image.
+- `proj4`: Projection definition (optional).
+- `wkt`: Well-known text definition (optional).
+- `epsg`: EPSG code of the image (optional).
 - `bbox`: The bounding boxes as a vector of GMTdataset.
 - `pixel_list`: A vector of vectors with the list of linear pixel indices for each component.
-- `centroid`: A vector of Float64 Tuples with the x,y coordinates of the centroids for each component.
+- `centroid`: A Float64 Matrix with the x,y coordinates of the centroids for each component.
 - `area`: A vector of Float64 with the areas of each component. This area is approximated by the mean lat for geog images.
 """
 bwconncomp(mat::BitMatrix; conn::Int=8) = bwconncomp(mat2img(collect(mat)); conn=conn)

src/spatial_funs.jl

@@ -178,13 +178,16 @@ 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`, `_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 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). If this attrib name is not found, we search for `CODE` which is
+  For that, name the keyword with a leading underscore, e.g. `_region`, `_nps`, `_area`, `_aspect`, `_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`,
+  `_aspect=0.5`. Areas are in square km when input is in geographic coordinates, otherwise square data units.
+  The aspect ratio passed to the `_aspect` option is width/height of the bounding box, not of the polygon itself.
+  Values <= 1, == 1 and >= 1 select polygons with aspect ratios less than, equal to and greater than the specified value.
+  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).
+  If this attrib name is not found, we search for `CODE` which is
   the one assigned by GMT when extracting polygons from the internal GMT coasts database. If none of these is found,
   it is users responsibility to provide a valid attribute name. The uniqueness is determined by selecting the polygon
   with the largest area.
@@ -210,7 +213,7 @@ Or `nothing` if the query results in an empty GMTdataset
 function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 	# This method does the work but it's not the one normally used by the public.
 	# It returns the indices of the selected segments.
-	(isempty(D[1].attrib)) && (@warn("This datset does not have an `attrib` field and is hence unusable here."); return Int[])
+	(isempty(D[1].attrib) && !is_in_kwargs(kw, [:_aspect])) && (@warn("This datset does not have an `attrib` field and is hence unusable here."); return Int[])
 	dounion = Int(1e9)		# Just a big number
 	invert = (find_in_kwargs(kw, [:invert :not :revert :reverse])[1] !== nothing)
 	if ((_att = find_in_kwargs(kw, [:att :attrib])[1]) !== nothing)		# For backward compat.
@@ -276,6 +279,11 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 		for k = 1:numel(D)  _tf[k] = areas_[k] >= area_  end
 		_tf
 	end
+	function clip_aspect(D, ratio, _tf)			# Clip by BBox aspect ratio
+		fun = (ratio >= 1) ? (>=) : (ratio == 1) ? (==) : <=
+		for k = 1:numel(D)  dd = diff(D[k].bbox); _tf[k] = fun(dd[1] / dd[3], ratio)  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.")
@@ -300,6 +308,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 .== "_aspect")) !== nothing && (ratio = parse(Float64, vals[ind]))
 	((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
 
@@ -312,6 +321,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] == "_aspect") tf = clip_aspect(D, ratio, tf)
 			elseif (atts[n] == "_unique") tf = clip_unique(D, areas, tf, attrib_name)
 			else                          tf = clip_np(D, nps, tf)
 			end
@@ -346,6 +356,7 @@ end
 See `? getbyattrib`
 """
 Base.:filter(D::Vector{<:GMTdataset}; kw...) = getbyattrib(D; kw...)
+Base.:filter(D::Vector{<:GMTdataset}, index::Bool; kw...) = getbyattrib(D, index; kw...)
 
 # ---------------------------------------------------------------------------------------------------
 """

src/utils.jl

@@ -111,6 +111,20 @@ function count_chars(str::AbstractString, c::Char=',')::Int
 	count(i->(i == c), str)
 end
 
+# ----------------------------------------------------------------------------------------------------------
+"""
+    ind = ind2bool(x::Vector{<:Integer} [, maxind::Int=0]) -> BitVector
+
+Convert a vector of indices `x` to a boolean vector `ind` where positions in `x` are true.
+
+If `maxind` is provided, the output boolean vector will have length `maxind`, otherwise its length will be `maximum(x)`.
+"""
+function ind2bool(indices::Vector{<:Integer}, maxind::Int=0)::BitVector
+	bool_vec = (maxind > 0) ? falses(maxind) : falses(maximum(indices))
+	bool_vec[indices] .= true
+	return bool_vec
+end
+
 # ----------------------------------------------------------------------------------------------------------
 """
     ind = uniqueind(x)