Merge pull request #1435 from GenericMappingTools/improve-polygonize

Change type of global POSTMAN. Improve polygonize to restrict size and order of output.

src/GMT.jl

@@ -86,7 +86,7 @@ const global noGrdCopy   = Vector{Bool}(undef, 1);noGrdCopy[1] = false		# If tru
 const global GMTCONF     = Vector{Bool}(undef, 1);GMTCONF[1] = false		# Flag if gmtset was used and must be 'unused' 
 const global FMT = ["png"]                         # The default plot format
 const global BOX_STR = [""]                        # Used in plotyy to know -R of first call
-const global POSTMAN = Dict{String,String}()       # To pass messages to functions (start with get_dataset) 
+const global POSTMAN = [Dict{String,String}()]     # To pass messages to functions (start with get_dataset) 
 const DEF_FIG_SIZE  = "15c/10c"                    # Default fig size for plot like programs. Approx 16/11
 const DEF_FIG_AXES_BAK     = " -Baf -BWSen"        # Default fig axes for plot like programs
 const DEF_FIG_AXES3_BAK    = " -Baf -Bza"          # 		"" but for 3D views

src/gdal_tools.jl

@@ -55,7 +55,7 @@ end
 Fill selected raster regions by interpolation from the edges.
 
 ### Parameters
-- `indata`:  Input data. It can be a file name, a GMTgrid or GMTimage object.
+- `data`:  Input data. It can be a file name, a GMTgrid or GMTimage object.
 - `nodata`: The nodata value that will be used to fill the regions. Otherwise use the `nodata` attribute of `indata`
    if it exists, or NaN if none of the above were set.
 - `kwargs`:
@@ -79,7 +79,7 @@ end
 Fill selected raster regions by interpolation from the edges.
 
 ### Parameters
-- `indata`:  Input data. The file name of a grid or image that can be read with `gmtread`.
+- `data`:  Input data. The file name of a grid or image that can be read with `gmtread`.
 - `nodata`: The nodata value that will be used to fill the regions. Otherwise use the `nodata` attribute of `indata`
    if it exists, or NaN if none of the above were set.
 - `kwargs`:
@@ -103,17 +103,33 @@ end
 This method, which uses the GDAL GDALPolygonize function, creates vector polygons for all connected regions
 of pixels in the raster sharing a common pixel/cell value. The input may be a grid or an image. This function can
 be rather slow as it picks lots of polygons in pixels with slightly different values at transitions between colors.
+Its natural use is to digitize masks images.
+
+### Parameters
+- `data`: Input data. It can be a GMTgrid or GMTimage object.
+- `kwargs`:
+  - `min, nmin, npixels or ncells`: The minimum number of cells/pixels for a polygon to be retained.
+    Default is 1. This can be set to filter out small polygons.
+  - `sort`: If true, will sort polygons by pixel count. Default is the order that GDAL decides internally.
 """
 function polygonize(data::GItype; gdataset=nothing, kwargs...)
 	d = GMT.KW(kwargs)
 	(gdataset === nothing) && (d[:gdataset] = true)
-	#(eltype(data) <: AbstractFloat)) && (d[:float] = true)	# To know which GDAL function to use.
+	if ((val = find_in_dict(d, [:min :nmin :npixels :ncells])[1]) !== nothing)
+		# Compute the cell area but ignoring eventual projection of if data is in geogs.
+		cell_area = Float64(val)::Float64 * data.inc[1] * data.inc[2];		s_area = string(cell_area)
+		isempty(GMT.POSTMAN[1]) ? (GMT.POSTMAN[1] = Dict("min_polygon_area" => s_area)) : GMT.POSTMAN[1]["min_polygon_area"] = s_area
+	end
 	o = helper_run_GDAL_fun(gdalpolygonize, data, "", String[], "", d...)
-	if (gdataset === nothing)		# Should be doing this for GDAL objects too but need to learn how to.
+	if (gdataset === nothing)						# Should be doing this for GDAL objects too but need to learn how to.
+		GMT.POSTMAN[1]["polygonize"] = "y"			# To inform gd2gmt() that it should check if last Di is the whole area.
+		(find_in_dict(d, [:sort])[1] !== nothing) && (GMT.POSTMAN[1]["sort_polygons"] = "y")
 		prj = getproj(data)
 		D = gd2gmt(o);		isa(D, Vector) ? (D[1].proj4 = prj) : (D.proj4 = prj)
+		delete!(GMT.POSTMAN[1], "min_polygon_area")		# In case it was set above
 		return D
 	end
+	delete!(GMT.POSTMAN[1], "min_polygon_area")
 	o
 end
 function polygonize(data::String; kwargs...)

src/gdal_utils.jl

@@ -254,7 +254,9 @@ function gd2gmt(dataset::Gdal.AbstractDataset)
 	# This method is for OGR formats only
 	(Gdal.OGRGetDriverByName(Gdal.shortname(getdriver(dataset))) == C_NULL) && return gd2gmt(dataset; pad=0)
 
-	D, ds = Vector{GMTdataset}(undef, Gdal.ngeom(dataset)), 1
+	min_area = (get(POSTMAN[1], "min_polygon_area", "") != "") ? parse(Float64, POSTMAN[1]["min_polygon_area"]) : 0.0
+	D, ds, get_area = Vector{GMTdataset}(undef, Gdal.ngeom(dataset)), 1, false
+	(get(POSTMAN[1], "sort_polygons", "") != "") && (polyg_area = zeros(length(D));		get_area = true)
 	proj = ""		# Fk local vars inside for 
 	for k = 1:Gdal.nlayer(dataset)
 		layer = getlayer(dataset, k-1)
@@ -267,8 +269,10 @@ function gd2gmt(dataset::Gdal.AbstractDataset)
 
 			for j = 0:Gdal.ngeom(feature)-1
 				geom = Gdal.getgeom(feature, j)
+				(min_area > 0.0 && geomarea(geom) < min_area) && continue
 				_D::GDtype = gd2gmt(geom, proj)
 				gt = Gdal.getgeomtype(geom)
+				get_area && (polyg_area[ds] = geomarea(geom))	# This area will NOT be what is expected if geom is known to be geog
 				# Maybe when there nlayers > 1 or other cases, starting allocated size is not enough
 				len_D = isa(_D, GMTdataset) ? 1 : length(_D)
 				(len_D + ds >= length(D)) && append!(D, Vector{GMTdataset}(undef, round(Int, 0.5 * length(D))))
@@ -293,6 +297,17 @@ function gd2gmt(dataset::Gdal.AbstractDataset)
 	D[1].colnames = startswith(proj, "+proj=longlat") ? ["lon","lat", ["z$i" for i=1:size(D[1].data,2)-2]...] :
 	                ["x","y", ["z$i" for i=1:size(D[1].data,2)-2]...]	
 	set_dsBB!(D)				# Compute and set the global BoundingBox for this dataset
+	if (get(POSTMAN[1], "polygonize", "") != "") && isapprox(D[1].ds_bbox, D[end].bbox)	# Last one is often an error (the whole area)
+		pop!(D)
+		delete!(GMT.POSTMAN[1], "polygonize")
+		ds -= 1
+	end
+	if (get(POSTMAN[1], "sort_polygons", "") != "")		# polygonize requested that the polygons go out in growing order
+		(polyg_area = deleteat!(polyg_area, ds:length(polyg_area)))
+		ind = sortperm(polyg_area, rev=true)
+		D = D[ind]
+		delete!(GMT.POSTMAN[1], "sort_polygons")
+	end
 	return (length(D) == 1) ? D[1] : D
 end
 

src/gmt_main.jl

@@ -609,10 +609,10 @@ function get_dataset(API::Ptr{Nothing}, object::Ptr{Nothing})::GDtype
 	D::GMT_DATASET = unsafe_load(convert(Ptr{GMT_DATASET}, object))
 
 	# This is for the particular case of the DCW countries that have a myriad of small segments and no Attributes
-	min_pts = (get(POSTMAN, "minpts", "") != "") ? parse(Int, POSTMAN["minpts"]) - 1 : 0
-	(min_pts > 0) && delete!(POSTMAN, "minpts")
-	DCWnames = (get(POSTMAN, "DCWnames", "") != "") ? true : false		# If DCW country names will turn into attribs
-	(DCWnames) && delete!(POSTMAN, "DCWnames")
+	min_pts = (get(POSTMAN[1], "minpts", "") != "") ? parse(Int, POSTMAN[1]["minpts"]) - 1 : 0
+	(min_pts > 0) && delete!(POSTMAN[1], "minpts")
+	DCWnames = (get(POSTMAN[1], "DCWnames", "") != "") ? true : false		# If DCW country names will turn into attribs
+	(DCWnames) && delete!(POSTMAN[1], "DCWnames")
 
 	seg_out = 0;
 	T::Vector{Ptr{GMT_DATATABLE}} = unsafe_wrap(Array, D.table, D.n_tables)

src/gmtreadwrite.jl

@@ -363,14 +363,14 @@ function file_has_time!(fname::String, D::GDtype, corder::Vector{Int}=Int[])
 end
 
 # ---------------------------------------------------------------------------------
-function guess_T_from_ext(fname::String, write::Bool=false)::String
+function guess_T_from_ext(fname::String; write::Bool=false, text_only::Bool=false)::String
 	# Guess the -T option from a couple of known extensions
 	fn, ext = splitext(fname)
 	if (ext == ".zip")			# Accept ogr zipped files, e.g., *.shp.zip
 		((out = guess_T_from_ext(fn)) == " -To") && return " -Toz"
 	end
 
-	_kml = (!write) ? "kml" : "*"	# This because on write we dont want to check for kml (let it be written as text)
+	_kml = (!write || !text_only) ? "kml" : "*"		# When it's text_only, we are writting an output gmt2kml
 
 	(length(ext) > 8 || occursin("?", ext)) && return (occursin("?", ext)) ? " -Tg" : "" # A SUBDATASET encoded fname?
 	ext = lowercase(ext[2:end])
@@ -380,7 +380,7 @@ function guess_T_from_ext(fname::String, write::Bool=false)::String
 	if     (findfirst(isequal(ext), ["grd", "nc", "nc=gd"])  !== nothing)  out = " -Tg";
 	elseif (findfirst(isequal(ext), ["dat", "txt", "csv"])   !== nothing)  out = " -Td";
 	elseif (findfirst(isequal(ext), ["jpg", "jpeg", "png", "bmp", "webp"]) 	!== nothing)  out = " -Ti";
-	elseif (findfirst(isequal(ext), ["arrow", "shp", _kml, "json", "feather", "geojson", "gmt", "gpkg", "gpx", "gml", "parquet"]) !== nothing)  out = " -To";
+	elseif (findfirst(isequal(ext), ["arrow", "shp", _kml, "kmz", "json", "feather", "geojson", "gmt", "gpkg", "gpx", "gml", "parquet"]) !== nothing)  out = " -To";
 	elseif (ext == "jp2") ressurectGDAL(); out = (findfirst("Type=UInt", gdalinfo(fname)) !== nothing) ? " -Ti" : " -Tg"
 	elseif (ext == "cpt")  out = " -Tc";
 	elseif (ext == "ps" || ext == "eps")  out = " -Tp";
@@ -524,8 +524,11 @@ function gmtwrite(fname::AbstractString, data; kwargs...)
 	(dbg_print_cmd(d, cmd) !== nothing) && return "gmtwrite " * fname * cmd
 
 	(opt_T == " -Tg" && isa(data, GMTgrid) && (data.scale != 1 || data.offset != 0)) && (fname *= "+s$(data.scale)+o$(data.offset)")
-	if (guess_T_from_ext(fname, true) == " -To")  gdalwrite(fname, data)		# Write OGR data
-	else                                          gmt("write " * fname * cmd, data)
+	text_only = isa(data, GMTdataset) && isempty(data.data)
+	if (guess_T_from_ext(fname, write=true, text_only=text_only) == " -To")
+		gdalwrite(fname, data)		# Write OGR data
+	else
+		gmt("write " * fname * cmd, data)
 	end
 	(opt_T == " -Ti") && transpcmap!(data, false)		# Reset original cmap (in case it was changed)
 	return nothing

src/pscoast.jl

@@ -94,10 +94,12 @@ function coast(cmd0::String=""; clip=nothing, first=true, kwargs...)
 
 	cmd = parse_E_coast(d, [:E, :DCW], "")		# Process first to avoid warning about "guess"
 	if (cmd != "")								# Check for a minimum of points that segments must have
-		((val = find_in_dict(d, [:minpts])[1]) !== nothing) && (POSTMAN["minpts"] = string(val)::String)
+		if ((val = find_in_dict(d, [:minpts])[1]) !== nothing)  POSTMAN[1]["minpts"] = string(val)::String
+		elseif (get(POSTMAN[1], "minpts", "") != "")            delete!(POSTMAN[1], "minpts")
+		end
 	end
 	cmd = add_opt(d, cmd, "M", [:M :dump])
-	contains(cmd, " -M") && (POSTMAN["DCWnames"] = "s")		# When dumping, we want to add the country name as attribute
+	contains(cmd, " -M") && (POSTMAN[1]["DCWnames"] = "s")		# When dumping, we want to add the country name as attribute
 	if (!occursin("-E+l", cmd) && !occursin("-E+L", cmd))
 		cmd, = parse_R(d, cmd, O)
 		if (!contains(cmd, " -M"))				# If Dump no -R & -B

test/test_gdal.jl

@@ -304,8 +304,8 @@ Gdal.GDALDestroyDriverManager()
 	fillnodata!(G);
 	@test GMT.isgeog(4326) == true
 
-	gmtwrite("9px.tif", mat2img(rand(UInt8, 3,3)))
-	polygonize("9px.tif");
+	gmtwrite("9px.tif", mat2img(round.(UInt8,round.(rand(6,6))*255)))
+	polygonize("9px.tif", nmin=2, sort=true);
 
 	gdaldrivers(:vec, out=true);
 end