diff --git a/src/gmt_main.jl b/src/gmt_main.jl
index 198d4eb74..ef709cafd 100644
--- a/src/gmt_main.jl
+++ b/src/gmt_main.jl
@@ -624,8 +624,11 @@ function get_dataset(API::Ptr{Nothing}, object::Ptr{Nothing})::GDtype
 	if (!isempty(POSTMAN[1]))
 		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")
+		if ((t = get(POSTMAN[1], "DCWnames", "")) != "")	# If DCW country names will turn into attribs
+			codelen = parse(Int, t)
+			delete!(POSTMAN[1], "DCWnames")
+		end
+		DCWnames = (t !== "") ? true : false
 		plusZzero = (get(POSTMAN[1], "plusZzero", "") != "") ? true : false		# To eventually add an extra column with 0's
 		(plusZzero) && delete!(POSTMAN[1], "plusZzero")
 	else
@@ -687,7 +690,7 @@ function get_dataset(API::Ptr{Nothing}, object::Ptr{Nothing})::GDtype
 					Darr[seg_out].header = hdrstr
 				else
 					(DCWnames && (ind = findfirst(" Segment", hdrstr)) !== nothing) &&
-						(Darr[seg_out].attrib["NAME"] = hdrstr[2:ind[1]-1])
+						(Darr[seg_out].attrib["CODE"] = hdrstr[4:4+codelen-1]; Darr[seg_out].attrib["NAME"] = hdrstr[7+codelen-2:ind[1]-1])
 				end
 			end
 			if (seg == 1)
diff --git a/src/pscoast.jl b/src/pscoast.jl
index 535b4a11b..3bb61c966 100644
--- a/src/pscoast.jl
+++ b/src/pscoast.jl
@@ -94,6 +94,7 @@ function _coast(cmd0::String, O::Bool, K::Bool, clip::String, d::Dict)
 	cmd = add_opt(d, "", "M", [:M :dump])
 	have_opt_M = (cmd != "")
 	cmd = parse_E_coast(d, [:E, :DCW], cmd)		# Process first to avoid warning about "guess"
+	twoORfour = have_opt_M && contains(cmd, "+z") && contains(cmd, '.') ? "4" : "2"		# To use in gmt_main to decide CODE attrib
 	if (cmd != "")								# Check for a minimum of points that segments must have
 		if ((val = find_in_dict(d, [:minpts])[1]) !== nothing)  POSTMAN[1]["minpts"] = string(val)::String
 		elseif (get(POSTMAN[1], "minpts", "") != "")            delete!(POSTMAN[1], "minpts")
@@ -104,7 +105,7 @@ function _coast(cmd0::String, O::Bool, K::Bool, clip::String, d::Dict)
 	toTrack::Union{String, GMTgrid} = ""
 	if (have_opt_M)
 		O = true
-		POSTMAN[1]["DCWnames"] = "s"		# When dumping, we want to add the country name as attribute
+		POSTMAN[1]["DCWnames"] = twoORfour		# When dumping, we want to add the country name as attribute
 		if ((val = find_in_dict(d, [:Z])[1]) !== nothing)
 			toTrack = (isa(val, GMTgrid) || (isa(val, String) && length(val) > 4)) ? val : ""
 			toTrack == "" && (POSTMAN[1]["plusZzero"] = "y")# If toTrack the extra column is added by the grdtrack call below
@@ -213,12 +214,13 @@ end
 function parse_E_coast(d::Dict, symbs::Vector{Symbol}, cmd::String)
 	(SHOW_KWARGS[1]) && return print_kwarg_opts(symbs, "NamedTuple | Tuple | Dict | String")
 	if ((val = find_in_dict(d, symbs, false)[1]) !== nothing)
-		if (isa(val, String) || isa(val, Symbol))	# Simple case, ex E="PT,+gblue" or E=:PT
+		if (isa(val, StrSymb))						# Simple cases, ex E="PT,+gblue" or E=:PT
 			t::String = string(" -E", val)
 			startswith(t, " -EWD") && (t = " -E=" * t[4:end])		# Let E="WD" work
 			(t == " -E") && (delete!(d, [:E, :DCW]); return cmd)	# This lets use E="" like earthregions may do
 			(!contains(cmd, " -M") && is_in_dict(d, [:R :region :limits :region_llur :limits_llur :limits_diag :region_diag]) === nothing) &&
 				(d[:R] = t[4:end])					# Must also see what to do for the other elseif branches
+			contains(cmd, " -M") && !contains(t, "+") && (t *= "+z")# If Dump always add country names as attribs
 			!contains(t, "+") && (t *= "+p0.5")		# If only code(s), append pen
 			cmd *= t
 		elseif (isa(val, NamedTuple) || isa(val, AbstractDict))
diff --git a/src/spatial_funs.jl b/src/spatial_funs.jl
index a0800ea7d..eccf467b1 100644
--- a/src/spatial_funs.jl
+++ b/src/spatial_funs.jl
@@ -184,7 +184,9 @@ NOTE: Instead of ``getbyattrib`` one case use instead ``filter`` (...,`index=fal
   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
+  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.
 
 - `invert, or reverse, or not`: If `true` return all segments that do NOT match the query condition(s).
@@ -245,6 +247,9 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 				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])
+					if (get(D[1].attrib, attrib_name, "") === "")
+						attrib_name = (get(D[1].attrib, "CODE", "") !== "") ? "CODE" : error("Attribute name '$attrib_name' or 'CODE' not found in dataset.")
+					end
 					vals[kk] = "0"	# If _unique the value doesn't matter but must be parseable to float
 				end
 				kk += 1
@@ -276,13 +281,19 @@ function getbyattrib(D::Vector{<:GMTdataset}, ind_::Bool; kw...)::Vector{Int}
 		((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.
+
+		if (numel(ind) == 1)					# The single group case is simpler but must be dealt separately.
+			 _tf[1] = false; _tf[argmax(areas_)] = true
+			 return _tf
+		end
+
 		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
+			_tf[n_start], _tf[this_ind] = false, true	# Turn off the first and turn on the largest of this group
 		end
 		_tf
 	end
diff --git a/test/test_B-GMTs.jl b/test/test_B-GMTs.jl
index 22e2ff233..58a52c27c 100644
--- a/test/test_B-GMTs.jl
+++ b/test/test_B-GMTs.jl
@@ -148,6 +148,7 @@
 	@test gmtselect([1 2], C="aa+d0", Vd=dbg2) == "gmtselect  -Caa+d0"
 	@test gmtselect([1 2], C=(pts=[1 2],dist=10), L=(line=[1 2;3 4], dist=10), Vd=dbg2) == "gmtselect  -C+d10 -L+d10"
 	cl = coast(dump=:true, res=:c, region=:global);
+	@test length(coast(M=true, DCW="WD+z")) > 0
 	cc = GMT.clipbyrect(cl, (-180,0,-90,90));
 
 	println("	GMTSET")