Merge pull request #1424 from GenericMappingTools/GI-empties

Need to cast the UInts floats in UInts GMTimage algebras

src/gdal_utils.jl

@@ -366,8 +366,8 @@ function gmt2gd(GI)
 	# Julia show methods. A further problem is that if we use 'GI = gmtread("file", layout="TRB")', GI's
 	# array still has a visible mem layout equal to a column major array. AS an example of this, see:
 	# gmtwrite("z.grd", mat2grid(reshape(1:20, 4,5))); Gr = gmtread("z.grd", layout="TRB"); Gc = gmtread("z.grd");
-	# Now both Gr.z and Gc.r are 4x5 matrices, but Gr can't be passed as is to 'gmt2gd(...)' because its
-	# memory layout is wrong when seen from the GDAL wrapper functions.  To disambiguate this, we are not
+	# Now both Gr.z and Gc.z are 4x5 matrices, but Gr can't be passed as is to 'gmt2gd(...)' because its
+	# memory layout is wrong when seen from the GDAL wrapper functions. To disambiguate this, we are not
 	# relying in 'size(GI)' but on the length of the x,y coordinates vectors, that are assumed to always be correct. 
 
 	width, height = (GI.layout != "" && GI.layout[2] == 'C') ? (size(GI,2), size(GI,1)) : (length(GI.x), length(GI.y)) .- GI.registration

src/gmt_main.jl

@@ -483,10 +483,10 @@ function get_image(API::Ptr{Nothing}, object)::GMTimage
 			o = (nz == 1) ? (ny, nx) : (ny, nx, nz)
 			isBRP = startswith(layout, "BRP")
 			(nz == 1 && isBRP) && (layout = "BRPa")	# For 1 layer "BRBa" and "BRPa" is actualy the same.
-			(!isBRP) && @warn("Only 'I' for Images.jl and 'BRP' MEM layouts are allowed.")
+			(!isBRP && nz > 1) && @warn("Only 'I' for Images.jl and 'BRP' MEM layouts are allowed.")
 
 			# OK, the above is not always true. Image cubes are not pixel interleaved. But how to detect them?
-			(I.type > 1 || (nz == 2 || nz > 4)) && (layout = layout[1:2] * "B")	# Stil leaves out cases of uint8 cubes.
+			(I.type > 1 || (nz == 2 || nz > 4)) && (layout = layout[1:2] * "B" * layout[4])	# Stil leaves out cases of uint8 cubes.
 		end
 		t = reshape(unsafe_wrap(Array, data, ny * nx * nz), o)	# Apparently the reshape() creates a copy
 	end
@@ -923,7 +923,8 @@ function image_init(API::Ptr{Nothing}, Img::GMTimage)::Ptr{GMT_IMAGE}
 		toRP_pad(Img, img_padded, n_rows, n_cols, pad)
 		already_converted = true
 	else
-		Ib.data = pointer(Img.image)
+		isa(Img.image, BitMatrix) && (copy_data = collect(Img.image))		# BitMatrx cannot be sent to C and we may need a copy bellow
+		Ib.data = (!isa(Img.image, BitMatrix)) ? pointer(Img.image) : pointer(copy_data)
 		mem_owned_by_gmt = (pad == 0) ? false : true
 	end
 
@@ -946,9 +947,9 @@ function image_init(API::Ptr{Nothing}, Img::GMTimage)::Ptr{GMT_IMAGE}
 	unsafe_store!(I, Ib)
 
 	if (!already_converted && !startswith(Img.layout, "BRP"))
-		img = (mem_owned_by_gmt) ? img_padded : copy(Img.image)
+		img = (mem_owned_by_gmt) ? img_padded : copy(Img.image)		# This copy is a waste when not Change_layout. Needs revisit.
 		(size(img,3) > 2) && GMT_Change_Layout(API, GMT_IS_IMAGE, "BRP", 0, I, img);	# Convert to BRP. Not 100% on the > 2 though.
-		Ib.data = pointer(img)
+		Ib.data = !isa(Img.image, BitMatrix) ? pointer(img) : pointer(copy_data)
 		unsafe_store!(I, Ib)
 	end
 
@@ -1414,6 +1415,7 @@ Shows information about the `D` GMTdataset (or vector of them).
 Runs ``show(stdout, "text/plain", any)`` which prints all elements of `any`. Good for printing the entire vector or matrix.
 """
 function info(GI::GItype, showdata::Bool=true; data=true, full=false, crs::Bool=false)
+	isempty(GI) && return println("Empty object")
 	crs && return print_crs(GI)
 	(data != 1) && (showdata = false)
 	isa(GI, GMTimage) ? println("A GMTimage object with $(size(GI,3)) bands of type $(eltype(GI))") :

src/gmt_types.jl

@@ -33,28 +33,28 @@ The fields of this struct are:
 - `hasnans::Int=0`:                     0 -> "don't know"; 1 -> confirmed, "have no NaNs"; 2 -> confirmed, "have NaNs"
 """
 Base.@kwdef mutable struct GMTgrid{T<:Number,N} <: AbstractArray{T,N}
-	proj4::String
-	wkt::String
-	epsg::Int
-	geog::Int
-	range::Union{Vector{Float64}, Vector{Any}}
-	inc::Union{Vector{Float64}, Vector{Any}}
-	registration::Int
-	nodata::Union{Float64, Float32}
-	title::String
-	remark::String
-	command::String
-	cpt::String
-	names::Vector{String}
-	x::Array{Float64,1}
-	y::Array{Float64,1}
-	v::Union{Vector{<:Real}, Vector{String}, Vector{<:TimeType}}
-	z::Array{T,N}
-	x_unit::String
-	y_unit::String
-	v_unit::String
-	z_unit::String
-	layout::String
+	proj4::String=""
+	wkt::String=""
+	epsg::Int=0
+	geog::Int=0
+	range::Union{Vector{Float64}, Vector{Any}}=Float64[]
+	inc::Union{Vector{Float64}, Vector{Any}}=Float64[]
+	registration::Int=0
+	nodata::Union{Float64, Float32}=0.0
+	title::String=""
+	remark::String=""
+	command::String=""
+	cpt::String=""
+	names::Vector{String}=String[]
+	x::Array{Float64,1}=Float64[]
+	y::Array{Float64,1}=Float64[]
+	v::Union{Vector{<:Real}, Vector{String}, Vector{<:TimeType}}=String[]
+	z::Array{T,N}=Array{Float64,2}(undef,0,0)
+	x_unit::String=""
+	y_unit::String=""
+	v_unit::String=""
+	z_unit::String=""
+	layout::String=""
 	scale::Union{Float64, Float32}=1f0
 	offset::Union{Float64, Float32}=0f0
 	pad::Int=0
@@ -106,28 +106,28 @@ The fields of this struct are:
 - `layout::String`:             A four character string describing the image memory layout
 - `pad::Int`:                   When != 0 means that the array is placed in a padded array of PAD rows/cols
 """
-mutable struct GMTimage{T<:Union{Unsigned, Bool}, N} <: AbstractArray{T,N}
-	proj4::String
-	wkt::String
-	epsg::Int
-	geog::Int
-	range::Vector{Float64}
-	inc::Vector{Float64}
-	registration::Int
-	nodata::Float32
-	color_interp::String
-	metadata::Vector{String}
-	names::Vector{String}
-	x::Vector{Float64}
-	y::Vector{Float64}
-	v::Vector{Float64}
-	image::Array{T,N}
-	colormap::Vector{Int32}
-	labels::Vector{String}		# Labels of a Categorical CPT
-	n_colors::Int
-	alpha::Matrix{UInt8}
-	layout::String
-	pad::Int
+Base.@kwdef mutable struct GMTimage{T<:Union{Unsigned, Bool, BitMatrix}, N} <: AbstractArray{T,N}
+	proj4::String=""
+	wkt::String=""
+	epsg::Int=0
+	geog::Int=0
+	range::Vector{Float64}=Float64[]
+	inc::Vector{Float64}=Float64[]
+	registration::Int=0
+	nodata::Float32=0f0
+	color_interp::String=""
+	metadata::Vector{String}=String[]
+	names::Vector{String}=String[]
+	x::Vector{Float64}=Float64[]
+	y::Vector{Float64}=Float64[]
+	v::Vector{Float64}=Float64[]
+	image::AbstractArray{T,N}=Array{UInt8,2}(undef,0,0)
+	colormap::Vector{Int32}=Int32[]
+	labels::Vector{String}=String[]		# Labels of a Categorical CPT
+	n_colors::Int=0
+	alpha::Matrix{UInt8}=Array{UInt8,2}(undef,0,0)
+	layout::String=""
+	pad::Int=0
 end
 Base.size(I::GMTimage) = size(I.image)
 Base.getindex(I::GMTimage{T,N}, inds::Vararg{Int,N}) where {T,N} = I.image[inds...]
@@ -216,19 +216,19 @@ The fields of this struct are:
 - `epsg::Int`:                    EPSG projection code (Optional)
 - `geom::Integer`:                Geometry type. One of the GDAL's enum (wkbPoint, wkbPolygon, etc...)
 """
-mutable struct GMTdataset{T<:Real, N} <: AbstractArray{T,N}
-	data::Array{T,N}
-	ds_bbox::Vector{Float64}
-	bbox::Vector{Float64}
-	attrib::DictSvS
-	colnames::Vector{String}
-	text::Vector{String}
-	header::String
-	comment::Vector{String}
-	proj4::String
-	wkt::String
-	epsg::Int
-	geom::Union{UInt32, Int}	# 0->Unknown, 1->Point, 2->Line, 3->Polygon, 4->MultiPoint, 5->MultiLine, 6->MultiPolyg
+Base.@kwdef mutable struct GMTdataset{T<:Real, N} <: AbstractArray{T,N}
+	data::Array{T,N}=Array{Float64,2}(undef,0,0)
+	ds_bbox::Vector{Float64}=Float64[]
+	bbox::Vector{Float64}=Float64[]
+	attrib::DictSvS=DictSvS()
+	colnames::Vector{String}=String[]
+	text::Vector{String}=String[]
+	header::String=""
+	comment::Vector{String}=String[]
+	proj4::String=""
+	wkt::String=""
+	epsg::Int=0
+	geom::Union{UInt32, Int}=0	# 0->Unknown, 1->Point, 2->Line, 3->Polygon, 4->MultiPoint, 5->MultiLine, 6->MultiPolyg
 end
 Base.size(D::GMTdataset) = size(D.data)
 Base.getindex(D::GMTdataset{T,N}, inds::Vararg{Int,N}) where {T,N} = D.data[inds...]
@@ -286,8 +286,8 @@ GMTdataset(data::Array{Float32,2}, text::String) =
 	GMTdataset(data, Float64[], Float64[], DictSvS(), String[], [text], "", String[], "", "", 0, 0)
 GMTdataset(data::Array{Float32,2}) =
 	GMTdataset(data, Float64[], Float64[], DictSvS(), String[], String[], "", String[], "", "", 0, 0)
-GMTdataset() =
-	GMTdataset(Array{Float64,2}(undef,0,0), Float64[], Float64[], DictSvS(), String[], String[], "", String[], "", "", 0, 0)
+#GMTdataset() =
+	#GMTdataset(Array{Float64,2}(undef,0,0), Float64[], Float64[], DictSvS(), String[], String[], "", String[], "", "", 0, 0)
 
 struct WrapperPluto fname::String end
 

src/grd_operations.jl

@@ -11,7 +11,7 @@ function Base.:+(G1::GMTgrid, G2::GMTgrid)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	hasnans = (G1.hasnans == 0 && G2.hasnans == 0) ? 0 : ((G1.hasnans + G2.hasnans) == 2) ? 1 : 2
 	(eltype(G1.z) <: AbstractFloat && eltype(G2.z) <: AbstractFloat) ? (z = G1.z .+ G2.z) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] + G2.z[k]  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) + Float32(G2.z[k])  end)
 	G3 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, hasnans)
 	setgrdminmax!(G3)		# Also take care of NaNs
@@ -26,7 +26,7 @@ function Base.:+(G1::GMTgrid, shift::Real)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	_shift = convert(eltype(G1.z), shift)
 	(eltype(G1.z) <: AbstractFloat) ? (z = G1.z .+ _shift) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] + _shift  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) + _shift  end)
 	G2 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	            z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, G1.hasnans)
 	G2.range[5:6] .+= shift
@@ -43,7 +43,7 @@ function Base.:-(G1::GMTgrid, G2::GMTgrid)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	hasnans = (G1.hasnans == 0 && G2.hasnans == 0) ? 0 : ((G1.hasnans + G2.hasnans) == 2) ? 1 : 2
 	(eltype(G1.z) <: AbstractFloat && eltype(G2.z) <: AbstractFloat) ? (z = G1.z .- G2.z) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] - G2.z[k]  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) - Float32(G2.z[k])  end)
 	G3 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, hasnans)
 	setgrdminmax!(G3)		# Also take care of NaNs
@@ -57,7 +57,7 @@ function Base.:-(G1::GMTgrid, shift::Real)
 	_shift = convert(eltype(G1.z), shift)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	(eltype(G1.z) <: AbstractFloat) ? (z = G1.z .- _shift) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] - _shift  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) - _shift  end)
 	G2 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, G1.hasnans)
 	G2.range[5:6] .-= shift
@@ -74,7 +74,7 @@ function Base.:*(G1::GMTgrid, G2::GMTgrid)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	hasnans = (G1.hasnans == 0 && G2.hasnans == 0) ? 0 : ((G1.hasnans + G2.hasnans) == 2) ? 1 : 2
 	(eltype(G1.z) <: AbstractFloat && eltype(G2.z) <: AbstractFloat) ? (z = G1.z .* G2.z) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] * G2.z[k]  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) * Float32(G2.z[k])  end)
 	G3 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, hasnans)
 	setgrdminmax!(G3)		# Also take care of NaNs
@@ -90,7 +90,7 @@ function Base.:*(G1::GMTgrid, scale::Real)
 	_scale = convert(eltype(G1.z), scale)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	(eltype(G1.z) <: AbstractFloat) ? (z = G1.z .* _scale) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] * _scale  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) * _scale  end)
 	G2 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, G1.hasnans)
 	G2.range[5:6] .*= scale
@@ -106,7 +106,7 @@ function Base.:^(G1::GMTgrid, scale::AbstractFloat)
 	_scale = convert(eltype(G1.z), scale)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	(eltype(G1.z) <: AbstractFloat) ? (z = G1.z .^ _scale) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] ^ _scale  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) ^ _scale  end)
 	G2 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	            z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, G1.hasnans)
 	G2.range[5:6] .^= scale
@@ -119,7 +119,7 @@ function Base.:/(G1::GMTgrid, G2::GMTgrid)
 	if (size(G1.z) != size(G2.z))  error("Grids have different sizes, so they cannot be divided.")  end
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	(eltype(G1.z) <: AbstractFloat && eltype(G2.z) <: AbstractFloat) ? (z = G1.z ./ G2.z) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] / G2.z[k]  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) / Float32(G2.z[k])  end)
 	hasnans = any(isnan, z) ? 2 : 1
 	G3 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, hasnans)
@@ -134,7 +134,7 @@ function Base.:/(G1::GMTgrid, scale::Real)
 	_scale = convert(eltype(G1.z), scale)
 	epsg, geog, range, inc, registration, nodata, x, y, v, pad = dup_G_meta(G1)
 	(eltype(G1.z) <: AbstractFloat) ? (z = G1.z ./ _scale) :
-		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = G1.z[k] / _scale  end)
+		(z = Matrix{Float32}(undef, size(G1.z)); @inbounds for k = 1:numel(G1.z)  z[k] = Float32(G1.z[k]) / _scale  end)
 	G2 = GMTgrid(G1.proj4, G1.wkt, epsg, geog, range, inc, registration, nodata, "", "", "", "", G1.names, x, y, v,
 	             z, G1.x_unit, G1.y_unit, G1.v_unit, G1.z_unit, G1.layout, 1f0, 0f0, pad, G1.hasnans)
 	G2.range[5:6] ./= scale
@@ -218,9 +218,22 @@ Base.:>=(I::GMTimage{T}, val::Real) where T <: Unsigned = mat2img(collect(I.imag
 """
 Subtract two boolean/uint8 mask images. It applies the logical `I1 && !I2` operation. Inherit header parameters from I1 image
 """
+#=
 Base.:-(I1::GMTimage{<:Bool},  I2::GMTimage{<:Bool})  = helper_bool_img(I1, collect(I1.image .& .!I2.image))
 Base.:-(I1::GMTimage{<:UInt8}, I2::GMTimage{<:UInt8}) =
 	helper_bool_img(I1, collect(reinterpret(Bool, I1.image) .& .!reinterpret(Bool, I2.image)))
+=#
+
+Base.:-(I1::GMTimage{<:Bool},  I2::GMTimage{<:Bool})  = helper_bool_img(I1, sub8_layout(I1, I2))
+Base.:-(I1::GMTimage{<:UInt8}, I2::GMTimage{<:UInt8}) = helper_bool_img(I1, sub8_layout(I1, I2))
+function sub8_layout(I1, I2)
+	I1.layout[2] == I2.layout[2] &&
+		return eltype(I1) <: Bool ? collect(I1.image .& .!I2.image) : collect(reinterpret(Bool, I1.image) .& .!reinterpret(Bool, I2.image))
+
+	z1 = (eltype(I1) <: Bool) ? I1.image  : reinterpret(Bool, I1.image)
+	z2 = (eltype(I2) <: Bool) ? I2.image' : reinterpret(Bool, I2.image)'
+	return collect(z1 .& .!z2)
+end
 
 # ---------------------------------------------------------------------------------------------------
 """
@@ -289,6 +302,15 @@ function Base.:!(I::GMTimage{<:UInt8})
 	return Io
 end
 
+"""
+    I = togglemask(I::Union{GMTimage{<:Bool, 2}, GMTimage{<:UInt8, 2}}) -> GMTimage
+
+Convert between UInt8 and Boolean representations of the mask images. A new object is returned but
+the underlying mask matrix is kept unchanged, that is, not copied, but reinterpred in the other type.
+"""
+togglemask(I::GMTimage{<:Bool, 2}) = return mat2img(reinterpret(UInt8, I.image), I)
+togglemask(I::GMTimage{<:UInt8, 2}) = return mat2img(reinterpret(Bool, I.image), I)
+
 # ---------------------------------------------------------------------------------------------------
 Base.:+(add::T, D1::GMTdataset) where T<:AbstractArray = Base.:+(D1::GMTdataset, add)
 Base.:+(add::Real, D1::GMTdataset) = Base.:+(D1::GMTdataset, [add;;])