When creating vectors of GMTdataset, assign them a data type as well (#…

…1685)

* Test/demo file for the whittaker smoothing function.

* Try reduce Vector{GMTdatasets} type ambiguities.

* More these

* Start giving type when creating empty Vector{GMTdataset{Float64,2}}

* Another bunch

* More

* Break up the plt_txt_attrib! function in two methods.

src/GMT.jl

@@ -199,7 +199,7 @@ include("gdal.jl")
 include("gdal_utils.jl")
 include("proj_utils.jl")
 using GMT.Gdal
-const global MatGDsGd = Union{Matrix{<:AbstractFloat}, GMTdataset, Vector{GMTdataset}, Gdal.AbstractDataset}
+const global MatGDsGd = Union{Matrix{<:AbstractFloat}, GMTdataset, Vector{<:GMTdataset}, Gdal.AbstractDataset}
 const global CURRENT_CPT = [GMTcpt()]		# To store the current palette
 
 include("gmt_main.jl")

src/gdal_extensions.jl

@@ -11,18 +11,18 @@ a `GMTdataset` or an `AbstractDataset`.
          type, in which case its referencing system is copied into `type`
 """
 function setproj!(tipo::AbstractArray, proj::String="")
-	(!isa(tipo, GMTgrid) && !isa(tipo, GMTimage) && !isa(tipo, GMTdataset) && !isa(tipo, Vector{GMTdataset})) &&
+	(!isa(tipo, GMTgrid) && !isa(tipo, GMTimage) && !isa(tipo, GMTdataset) && !isa(tipo, Vector{<:GMTdataset})) &&
 		error("Wrong data type for this function. Must be a grid, image or dataset")
 	(proj == "") && error("the projection string cannot obviously be empty")
 	isproj4 = (startswith(proj, "+proj") !== nothing)
-	obj = (isa(tipo, Vector{GMTdataset})) ? tipo[1] : tipo
+	obj = (isa(tipo, Vector{<:GMTdataset})) ? tipo[1] : tipo
 	(isproj4) ? (obj.proj4 = proj) : (obj.wkt = proj)
 	return nothing
 end
 function setproj!(tipo::AbstractArray, ref)
-	(!isa(ref, GMTgrid) && !isa(ref, GMTimage) && !isa(ref, GMTdataset) && !isa(ref, Vector{GMTdataset})) &&
+	(!isa(ref, GMTgrid) && !isa(ref, GMTimage) && !isa(ref, GMTdataset) && !isa(ref, Vector{<:GMTdataset})) &&
 		error("Wrong REFERENCE data type for this function. Must be a grid, image or dataset")
-	obj = (isa(ref, Vector{GMTdataset})) ? ref[1] : ref
+	obj = (isa(ref, Vector{<:GMTdataset})) ? ref[1] : ref
 	((prj = obj.proj4) == "") && (prj = obj.wkt)
 	(prj == "") && error("The REFERENCE type is not referenced with either PROJ4 or WKT string")
 	setproj!(tipo, prj)

src/gdal_utils.jl

@@ -219,7 +219,7 @@ function gd2gmt(geom::Gdal.AbstractGeometry, proj::String="")::Union{GMTdataset,
 		geom = Gdal.getgeom(geom,0)
 	elseif (gmtype == wkbMultiPolygon || gmtype == wkbMultiLineString || gmtype == Gdal.wkbGeometryCollection || gmtype == Gdal.wkbGeometryCollection25D)
 		n_pts = Gdal.ngeom(geom)
-		D = Vector{GMTdataset}(undef, n_pts)
+		D = Vector{GMTdataset{Float64, 2}}(undef, n_pts)
 		for k = 1:n_pts  D[k] = gd2gmt(Gdal.getgeom(geom,k-1), "")::GMTdataset  end
 		(proj != "") && (D[1].proj4 = proj)
 		set_dsBB!(D)				# Compute and set the BoundingBox's for this dataset
@@ -262,7 +262,7 @@ function gd2gmt(dataset::Gdal.AbstractDataset)
 	min_area = (get(POSTMAN[1], "min_polygon_area", "") != "") ? parse(Float64, POSTMAN[1]["min_polygon_area"]) : 0.0
 	max_area = (get(POSTMAN[1], "max_polygon_area", "") != "") ? parse(Float64, POSTMAN[1]["max_polygon_area"]) : 0.0
 	p_isgeog = (get(POSTMAN[1], "polygon_isgeog", "") != "") ? true : false
-	D, ds, get_area = Vector{GMTdataset}(undef, Gdal.ngeom(dataset)), 1, false
+	D, ds, get_area = Vector{GMTdataset{Float64, 2}}(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)
@@ -283,7 +283,7 @@ function gd2gmt(dataset::Gdal.AbstractDataset)
 				get_area && (polyg_area[ds] = geomarea(geom))	# This area will NOT be what is expected if geom is known to be geog
 				# Maybe when 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))))
+				(len_D + ds >= length(D)) && append!(D, Vector{GMTdataset{Float64, 2}}(undef, round(Int, 0.5 * length(D))))
 				if isa(_D, GMTdataset)
 					D[ds] = _D
 					D[ds].geom = gt
@@ -354,7 +354,7 @@ end
 # would not be able to extract the data from the 'dataset'
 function helper_read_XLSCSV(dataset::Gdal.AbstractDataset)::GDtype
 	n_layers = Gdal.nlayer(dataset)
-	D = Vector{GMTdataset}(undef, n_layers)
+	D = Vector{GMTdataset{Float64, 2}}(undef, n_layers)
 	local inds_r, inds_s					# n^inf Fck. If I don't do this here errors are completly non-sensical.
 	for n_layer = 1:n_layers
 		layer = getlayer(dataset, n_layer-1)
@@ -1047,7 +1047,7 @@ Return a vector of GMTdatasets from a vector of WKT strings.
 """
 function readgeom(wkt::Vector{String})
 	n_ds = length(wkt)
-	D::Vector{GMTdataset} = Vector{GMTdataset}(undef, n_ds)
+	D = Vector{GMTdataset}(undef, n_ds)
 	for k = 1:n_ds
 		D[k] = gd2gmt(fromWKT(wkt[k]))
 	end

src/gmt_main.jl

@@ -1371,7 +1371,7 @@ function ogr2GMTdataset(in::Ptr{OGR_FEATURES}, drop_islands=false)::Union{GMTdat
 		(n_islands > 0) && println("\tThis file has islands (holes in polygons).\n\tUse `gmtread(..., no_islands=true)` to ignore them.")
 	end
 
-	D::Vector{GMTdataset} = Vector{GMTdataset}(undef, n_total_segments)
+	D = Vector{GMTdataset{Float64, 2}}(undef, n_total_segments)
 
 	n = 1
 	attrib = DictSvS();	# For the case there are no attribs at all.
@@ -1702,7 +1702,7 @@ Base.:show(io::IO, mime::MIME"image/png", wp::WrapperPluto) = write(io, read(wp.
 # ---------- For Tables -----------------------------------------------------------------------------
 Base.:names(D::GDtype) = isa(D, Vector) ? D[1].colnames : D.colnames
 #Base.:names(D::GMTdataset) = D.colnames
-#Base.:names(D::Vector{GMTdataset}) = D[1].colnames
+#Base.:names(D::Vector{<:GMTdataset}) = D[1].colnames
 
 # ---------- For fck stop printing UInts in hexadecinal ---------------------------------------------
 #Base.show(io::IO, x::T) where {T<:Union{UInt, UInt128, UInt64, UInt32, UInt16, UInt8}} = Base.print(io, x)

src/gmt_types.jl

@@ -342,7 +342,7 @@ Base.@kwdef struct GMTtypes
 	grd::GMTgrid = GMTgrid()
 	img::GMTimage = GMTimage()
 	ds::GMTdataset = GMTdataset()
-	dsv::Vector{GMTdataset} = [GMTdataset()]
+	dsv::Vector{GMTdataset::Vector{GMTdataset{Float64, 2}}} = [GMTdataset()]
 	cpt::GMTcpt = GMTcpt()
 	ps::GMTps = GMTps()
 	function GMTtypes(stored, grd, img, ds, dsv, cpt, ps)

src/gmtbegin.jl

@@ -473,7 +473,7 @@ function zoom2inset(d, center)
 	else
 		bak = CTRL.limits[7:end]						# Backup these because round_wesn will change them
 		nDS = (isa(data, Vector)) ? length(data) : 1
-		D = Vector{GMTdataset}(undef, nDS)
+		D = Vector{GMTdataset{Float64, 2}}(undef, nDS)
 		zoom_lims::Vector{Float64} = [Inf, -Inf, Inf, -Inf]
 		for k = 1:nDS
 			if (nDS == 1)

src/gmtselect.jl

@@ -108,7 +108,7 @@ function clipbyrect(D::Vector{<:GMTdataset}, limits; revert::Bool=false)
 		cross[k] = !in[k] && !out
 	end
 	revert && (in .= .!in)
-	Dclipped = Vector{GMTdataset}(undef, sum(in)+sum(cross))
+	Dclipped = Vector{GMTdataset{Float64,2}}(undef, sum(in)+sum(cross))
 	m = 0
 	for k = 1:numel(D)
 		in[k] && (Dclipped[m+=1] = deepcopy(D[k]))

src/grdinterpolate.jl

@@ -175,7 +175,7 @@ function grdinterp_local_opt_S(arg1::GItype, pts::Vector{<:GMTdataset}, no_coord
 	# This whole thing needs further documentation and testing.
 	n_layers = size(arg1,3)
 	startcol = no_coords ? 0 : 2
-	D = Vector{GMTdataset}(undef, length(pts))
+	D = Vector{GMTdataset{Float64,2}}(undef, length(pts))
 	DT = no_coords ? eltype(arg1) : eltype(pts[1])		# When we have coordinates, their type dominates.
 	(DT <: Integer) && (DT = eltype(pts[1]))			# We don't want conversions to integers that are a pain. Specially if UInt8.
 	for n = 1:length(pts)					# Loop over number of segments and initialize the output dataset vector

src/laszip/lazread.jl

@@ -2,7 +2,7 @@ Base.@kwdef struct lasout_types
 	stored::String = ""
 	grd::GMTgrid = GMTgrid()
 	ds::GMTdataset = GMTdataset()
-	dsv::Vector{GMTdataset} = [GMTdataset()]
+	dsv::Vector{GMTdataset{Float64}} = [GMTdataset()]
 	function lasout_types(stored, grd, ds, dsv)
 		stored = !isempty(grd) ? "grd" : (!isempty(ds) ? "ds" : (!isempty(dsv) ? "dsv" : ""))
 		new(stored, grd, ds, dsv)

src/lowess.jl

@@ -66,10 +66,7 @@ function lowess(x::AbstractVector{T}, y::AbstractVector{T}; span::T=2/3, nsteps:
 	c1::T = 0.0
 	r::T = 0.0
 
-	if (n < 2)
-		ys[1] = y[1]
-		return ys
-	end
+	(n < 2) && (ys[1] = y[1]; return ys)
 
 	ns = max(min(floor(Int, span * n), n), 2)  # at least two, at most n points
 	for iter = 1:(nsteps + 1)  # robustness iterations
@@ -93,9 +90,7 @@ function lowess(x::AbstractVector{T}, y::AbstractVector{T}; span::T=2/3, nsteps:
 			lowest(x, y, n, x[i + 1], ys, i, nleft, nright, res, (iter > 1), rw, ok)
 
 			# fitted value at x[i + 1]
-			if (ok == 0)
-				ys[i + 1] = y[i + 1]
-			end
+			(ok == 0) && (ys[i + 1] = y[i + 1])
 
 			# all weights zero - copy over value (all rw==0)
 			if (last < i - 1)   # skipped points -- interpolate

src/pca.jl

@@ -227,7 +227,7 @@ function kmeans(X::Union{GMTdataset, Matrix{<:Real}}, k=3; seeds=nothing, maxite
 
 	(raw) && return idx, centers, [classes [sum(idx .== classes[n]) for n = 1:numel(classes)]]
 
-	Dv = Vector{GMTdataset}(undef, length(classes))
+	Dv = Vector{GMTdataset{Float64,2}}(undef, length(classes))
 	colnames = isa(X, GMTdataset) ? X.colnames : String[]
 	has_text = isa(X, GMTdataset) && !isempty(X.text)
 	for n = 1:numel(classes)

src/pcolor.jl

@@ -104,7 +104,7 @@ function pcolor(X_::VMr, Y_::VMr, C::Union{Nothing, AbstractMatrix{<:Real}}, fir
 		[Y[k] -= yinc2 for k = 1:numel(Y)];	append!(Y, Y[end]+yinc)
 	end
 
-	D::Vector{GMTdataset}, k = Vector{GMTdataset}(undef, length(C)), 0
+	D, k = Vector{GMTdataset{Float64,2}}(undef, length(C)), 0
 	if (isvector(X))
 		for col = 1:length(X)-1, row = 1:length(Y)-1	# Gdal.wkbPolygon = 3
 			if (k == 0) 
@@ -283,15 +283,15 @@ function boxes(X::VMr, Y::VMr, d::Dict{Symbol,Any})
 	end
 
 	n_tiles = isvector(X) ? (length(X) - 1)*(length(Y) - 1) : (size(X,1) - 1)*(size(X,2) - 1)
-	D::Vector{GMTdataset}, k = Vector{GMTdataset}(undef, n_tiles), 0
+	D, k = Vector{GMTdataset{promote_type(eltype(X), eltype(Y)), 2}}(undef, n_tiles), 0
 	if (isvector(X))
 		for col = 1:length(X)-1, row = 1:length(Y)-1
 			(isautomask && Gmask.z[row, col] != mask_true) && continue
 			if (k == 0) D[k+=1] = mat2ds([X[col] Y[row]; X[col] Y[row+1]; X[col+1] Y[row+1]; X[col+1] Y[row]; X[col] Y[row]]; geom=3, d...)
 			else        D[k+=1] = mat2ds([X[col] Y[row]; X[col] Y[row+1]; X[col+1] Y[row+1]; X[col+1] Y[row]; X[col] Y[row]]; geom=3)
 			end
 		end
-		k == 0 && return GMTdataset[]
+		k == 0 && return GMTdataset{eltype(D.data),2}[]
 		(isautomask && k != (length(X)-1)*(length(Y)-1)) && deleteat!(D, k+1:n_tiles)	# Remove the unused D's
 		D[1].ds_bbox = [X[1], X[end], Y[1], Y[end]]
 		isautomask && set_dsBB!(D, false) 

src/plot.jl

@@ -109,7 +109,7 @@ function plot(arg1; first=true, kw...)
 		gidx, gnames = get_group_indices(d, arg1)
 		cycle_colors = (numel(gidx) <= 7) ? matlab_cycle_colors : simple_distinct	# Will blow if > 20
 		if (!isempty(gidx))
-			Dv = Vector{GMTdataset}(undef, length(gidx))
+			Dv = Vector{GMTdataset{Float64,2}}(undef, length(gidx))
 			for k = 1:numel(gidx)
 				Dv[k] = linearfitxy(mat2ds(arg1, (gidx[k], :)))
 				Dv[k].header = "-G"*cycle_colors[k] * " -W"*cycle_colors[k]
@@ -738,7 +738,7 @@ function fill_between(arg1, arg2=nothing; first=true, kwargs...)
 		end
 		ind = find_the_pos(view(D1, :, 1), view(int, :, 1))		# Indices of the points before the intersections
 		n_crossings = size(ind,1)
-		Dsd = Vector{GMTdataset}(undef, n_crossings+1)
+		Dsd = Vector{GMTdataset{Float64,2}}(undef, n_crossings+1)
 		ff = (D1[1,2] < D1[1,3]) ? 2 : 1
 		Dsd[1] = mat2ds([D1[1:ind[1], [1,2]]; int[1:1,1:2]; D1[ind[1]:-1:1, [1,3]]], fill=fill_colors[ff])
 		for k = 2:n_crossings
@@ -760,7 +760,7 @@ function fill_between(arg1, arg2=nothing; first=true, kwargs...)
 		ind1 = find_the_pos(view(D1, :, 1), view(int, :, 1))		# Indices of the points before the intersections at line 1
 		ind2 = find_the_pos(view(D2, :, 1), view(int, :, 1))		# Indices of the points before the intersections at line 2
 		n_crossings = size(int,1)
-		Dsd = Vector{GMTdataset}(undef, n_crossings+1)
+		Dsd = Vector{GMTdataset{Float64,2}}(undef, n_crossings+1)
 		ff = (D1[1,2] < D2[1,2]) ? 2 : 1
 		Dsd[1] = mat2ds([D1[1:ind1[1], [1,2]]; int[1:1,1:2]; D2[ind2[1]:-1:1, [1,2]]], fill=fill_colors[ff])
 		for k = 2:n_crossings
@@ -1531,7 +1531,7 @@ function helper_hvband(mat::Matrix{<:Real}, tipo="v"; width=false, height=false,
 	else              ind_w, ind_b, ind_t, thick, bB = 9:10, 7, 9, height != 0, "B"
 	end
 
-	D::Vector{GMTdataset} = Vector{GMTdataset}(undef, n_ds)
+	D = Vector{GMTdataset{Float64,2}}(undef, n_ds)
 	for k = 1:n_ds
 		w = (thick) ? mat[k,2] : (percent != 0) ? mat[k,2]*diff(CTRL.limits[ind_w]) : mat[k,2]-mat[k,1]	# bar width
 		i = rem(k, length(colors)); (i == 0) && (i = length(colors))

src/proj_utils.jl

@@ -312,7 +312,7 @@ buffergeo(D::GMTdataset; width=0, unit=:m, np=120, flatstart=false, flatend=fals
 	buffergeo(D.data; width=width, unit=unit, np=np, flatstart=flatstart, flatend=flatend, proj=D.proj4, epsg=epsg, tol=tol)
 
 function buffergeo(D::Vector{<:GMTdataset}; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=-1.0)
-	_D = Vector{GMTdataset}(undef, length(D))
+	_D = Vector{GMTdataset{Float64,2}}(undef, length(D))
 	for k = 1:length(D)
 		_D[k] = buffergeo(D[k].data; width=width, unit=unit, np=np, flatstart=flatstart, flatend=flatend,
 		                  proj=D[1].proj4, epsg=epsg, tol=tol)
@@ -432,7 +432,7 @@ geodesic(ds::Gdal.AbstractDataset; step=0.0, unit=:m, np=0, proj::String="", eps
 
 function geodesic(D::Vector{<:GMTdataset}; step=0.0, unit=:m, np=0, proj::String="", epsg::Integer=0,
 	              longest::Bool=false)
-	_D = Vector{GMTdataset}(undef, length(D))
+	_D = Vector{GMTdataset{Float64,2}}(undef, length(D))
 	for k = 1:length(D)
 		_D[k] = mat2ds(geodesic(D[k].data; step=step, unit=unit, np=np, proj = (proj == "") ? D[k].proj4 : proj,
 		               epsg=epsg, longest=longest))
@@ -613,7 +613,7 @@ function helper_gdirect(projPJ_ptr, lonlat, azim, dist, proj_string, isgeog, dat
 			(length(dist) != n_lines) && (proj_destroy(projPJ_ptr);error("Number of distance vectors MUST be equal to number of azimuths"))
 			Vdist = dist
 		end
-		D = Vector{GMTdataset}(undef, n_lines)
+		D = Vector{GMTdataset{Float64,2}}(undef, n_lines)
 		for nl = 1:n_lines
 			n_pts = length(Vdist[nl])					# Number of points in this line
 			dest = Array{Float64}(undef, n_pts, 3)		# Azimuth goes into the D too

src/psxy.jl

@@ -404,7 +404,8 @@ function parse_plot_callers(d, gmt_proggy, caller, is3D, O, arg1)
 end
 
 # ---------------------------------------------------------------------------------------------------
-function plt_txt_attrib!(D::GDtype, d::Dict, _cmd::Vector{String})
+plt_txt_attrib!(D::GMTdataset, d::Dict, _cmd::Vector{String}) = plt_txt_attrib!([D], d, _cmd)
+function plt_txt_attrib!(D::Vector{<:GMTdataset}, d::Dict, _cmd::Vector{String})
 	# Plot TEXT attributed labels and serve as function barrier agains to f Anys
 	((val = find_in_dict(d, [:labels])[1]) === nothing) && return nothing
 
@@ -1725,7 +1726,7 @@ because we may want to see the inside of a top surface.
 and `x` and `y` are in degrees, case in which an automatic projection takes place, or in meters. We need this
 if we want that the normal compuations makes sense.
 """
-function sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfact=1.0)::Vector{GMTdataset}
+function sort_visible_triangles(Dv::Vector{<:GMTdataset}; del_hidden=false, zfact=1.0)
 	azim, elev = parse.(Float64, split(CURRENT_VIEW[1][4:end], '/'))
 	sin_az, cos_az, sin_el = sind(azim), cosd(azim), sind(elev)
 	prj, wkt, epsg = Dv[1].proj4, Dv[1].wkt, Dv[1].epsg
@@ -1818,14 +1819,14 @@ or, to plot them
 viz(FV, replicate=(centers=rand(10,3)*10, scales=0.1))
 ```
 """
-function replicant(FV::GMTfv; kwargs...)::Vector{GMTdataset}		# For direct calls to replicat()
+function replicant(FV::GMTfv; kwargs...)		# For direct calls to replicat()
 	d = KW(kwargs)
 	(is_in_dict(d, [:p :view :perspective]) === nothing) && (CURRENT_VIEW[1] = " -p217.5/30")
 	replicant(FV, d)
 end
 
 # ---------------------------------------------------------------------------------------------------
-function replicant(FV::GMTfv, d::Dict{Symbol, Any})::Vector{GMTdataset}
+function replicant(FV::GMTfv, d::Dict{Symbol, Any})
 	(val = find_in_dict(d, [:replicate])[1]) === nothing && error("Can't replicate without the 'replicate' option")
 
 	cpt::GMTcpt = GMTcpt()
@@ -1867,16 +1868,16 @@ function replicant(FV::GMTfv, d::Dict{Symbol, Any})::Vector{GMTdataset}
 end
 
 # ---------------------------------------------------------------------------------------------------
-function replicant_worker(FV::GMTfv, xyz, azim, elev, cval, cpt, scales)::Vector{GMTdataset}
+function replicant_worker(FV::GMTfv, xyz, azim, elev, cval, cpt, scales)
 	# This guy is the one who does the replicant work
 
 	FV, normals = sort_visible_faces(FV, azim, elev)	# Return a modified FV containing info about the sorted visible faces.
 
 	n_faces_tot = sum(size.(FV.faces_view, 1))			# Total number of segments or faces (polygons)
 
 	# ---------- First we convert the FV into a vector of GMTdataset
-	D1 = Vector{GMTdataset}(undef, n_faces_tot)
 	t = zeros(eltype(FV.verts), maximum(size.(FV.faces_view,2)), 3)	# The maximum number of vertices in any polygon of all geometries
+	D1 = Vector{GMTdataset{eltype(t),2}}(undef, n_faces_tot)
 	count_face = 0
 
 	for geom = 1:numel(FV.faces_view)					# If we have more than one type of geometries (e.g. triangles and quadrangles)
@@ -1894,7 +1895,7 @@ function replicant_worker(FV::GMTfv, xyz, azim, elev, cval, cpt, scales)::Vector
 	P::Ptr{GMT_PALETTE} = palette_init(G_API[1], cpt)	# A pointer to a GMT CPT
 	rgb = [0.0, 0.0, 0.0, 0.0]
 	cor = [0.0, 0.0, 0.0]
-	D2 = Vector{GMTdataset}(undef, size(xyz, 1) * n_faces_tot)
+	D2 = Vector{GMTdataset{promote_type(eltype(xyz), eltype(scales), eltype(t)),2}}(undef, size(xyz, 1) * n_faces_tot)
 
 	# ---------- Now we do the replication
 	for k = 1:size(xyz, 1)								# Loop over number of positions. For each of these we have a new body

src/spatial_funs.jl

@@ -444,7 +444,7 @@ A GMTdatset if only one polygon was passed or a Vector{GMTaset} otherwise.
 randinpolygon(mat::Matrix{<:AbstractFloat}; density=0.1, np::Int=0) = randinpolygon(mat2ds(mat); density=density, np=np)
 function randinpolygon(Din::GDtype; density=0.1, np::Int=0)
 
-	D::Vector{GMTdataset} = Vector{GMTdataset}(undef, isa(Din, Vector) ? length(Din) : 1)
+	D = Vector{GMTdataset{Float64, 2}}(undef, isa(Din, Vector) ? length(Din) : 1)
 
 	function get_the_points(_D, dx, dy, density, np, isgeo)
 		# Points are returned in same data types as those of _D