Let gmtread and write deal with lasziped files as well.

src/GMT.jl

@@ -331,7 +331,7 @@ include("zscale.jl")
 include("drawing.jl")
 include("get_enums.jl")
 
-include("laszip/laszip.jl")
+include("laszip/Laszip.jl")
 using GMT.Laszip
 
 @setup_workload let

src/gmtreadwrite.jl

@@ -105,7 +105,7 @@ function gmtread(_fname::String; kwargs...)
 			if ((val1 = find_in_dict(d, [:layer :layers :band :bands])[1]) !== nothing)
 				if (isa(val1, Real))               gdopts = string(" -b ", val1)
 				elseif (isa(val1, AbstractArray))  gdopts = join([string(" -b ", val1[i]) for i in 1:numel(val1)])
-			end
+				end
 			end
 		else
 			fname *= "?" * varname
@@ -168,6 +168,9 @@ function gmtread(_fname::String; kwargs...)
 			fname, opt_T = "/vsizip/" * fname, " -To"
 		elseif (opt_T == "obj")						# Means we got a .obj file. Read it and leave
 			return read_obj(fname)
+		elseif (opt_T == "las")						# Means we got a .laz or .las file. Read it and leave
+			o_las = las2dat(fname; kwargs...)
+			return getproperty(o_las, Symbol(o_las.stored))
 		end
 	else
 		opt_T = opt_T[1:4]      					# Remove whatever was given as argument to type kwarg
@@ -422,6 +425,7 @@ function guess_T_from_ext(fname::String; write::Bool=false, text_only::Bool=fals
 	fn, ext = splitext(fname)
 	ext = lowercase(ext[2:end])
 	(ext == "obj") && return "obj"	# To be read by read_obj() internal function.
+	(ext == "laz" || ext == "las") && return "las"	# To be read by las2dat()
 	if (ext == "zip")				# Accept ogr zipped files, e.g., *.shp.zip
 		((out = guess_T_from_ext(fn)) == " -To") && return " -Toz"
 	end
@@ -517,10 +521,7 @@ function gmtwrite(fname::AbstractString, data; kwargs...)
 	(fname == "") && error("Output file name cannot be empty.")
 
 	if (isa(data, GMTgrid))
-		#opt_T = " -Tg"
-		#fname *= parse_grd_format(d)			# If we have format requests
-		#cmd, = parse_f(d, cmd)
-		#CTRL.proj_linear[1] = true				# To force pad=0 and julia memory (no dup)
+		(endswith(fname, ".laz") || endswith(fname, ".LAZ")) && return dat2las(fname, data; kwargs...)		# Lasz
 
 		# GMT doesn't write correct CF nc grids that are referenced but non-geographic. So, use GDAL in those cases
 		fmt = parse_grd_format(d)				# See if we have format requests
@@ -543,6 +544,7 @@ function gmtwrite(fname::AbstractString, data; kwargs...)
 		transpcmap!(data, true)
 	elseif (isa(data, GDtype))
 		isa(data, Vector) && (endswith(fname, ".stl") || endswith(fname, ".STL")) && return write_stl(fname, data; kwargs...)	# STL
+		(endswith(fname, ".laz") || endswith(fname, ".LAZ")) && return dat2las(fname, data; kwargs...)		# Lasz
 		opt_T = " -Td"
 		cmd, = parse_bo(d, cmd)					# Write to binary file
 	elseif (isa(data, GMTcpt))
@@ -560,6 +562,7 @@ function gmtwrite(fname::AbstractString, data; kwargs...)
 			opt_T = " -Ti"
 		end
 	elseif (isa(data, AbstractArray))
+		(endswith(fname, ".laz") || endswith(fname, ".LAZ")) && return dat2las(fname, data; kwargs...)		# Lasz
 		fmt = parse_grd_format(d)				# See if we have format requests
 		if (fmt == "")							# If no format, write a dataset
 			opt_T = " -Td"

src/laszip/Laszip.jl

@@ -1,6 +1,6 @@
 module Laszip
 
-using Printf, Dates, LASzip_jll
+using GMT, Printf, Dates, LASzip_jll
 
 export
 	xyz2laz, laz2xyz, las2dat, dat2las

src/laszip/dat2las.jl

@@ -1,30 +1,39 @@
 """
-Write XYZ data to a LIDAR laz (laszip compressed) or las format file. Usage:
+    dat2las(FileName::AbstractString, xyz; grd_hdr=[], scaleX=nothing, scaleY=nothing, scaleZ=nothing, offX=nothing, offY=nothing, offZ=nothing)
 
-	dat2las(FileName::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], scaleZ=[], offX=[], offY=[], offZ=[])
+Write XYZ data to a LIDAR laz (laszip compressed) or las format file.
 
 	Where:
 		"FileName" Name of the output LIDAR file
 		xyz  A Mx3 array with the point coordinates
 
-Example. To write the x,y,z data to file "lixo.laz" do:
+### Example
 
+To write the x,y,z data to file "lixo.laz" do:
+
+```julia
 	dat2las("lixo.laz", xyz)
+```
 """
-function dat2las(fname::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], scaleZ=[], offX=[], offY=[], offZ=[])
+function dat2las(fname::AbstractString, G::GMTgrid; scaleX=nothing, scaleY=nothing, scaleZ=nothing, offX=nothing, offY=nothing, offZ=nothing)
+	if     (startswith(G.layout, "BC"))  z = vec(G.z[:])
+	elseif (startswith(G.layout, "TR"))  z = grd2xyz(G, Z="TLf")
+	end
+	hdr = [G.range[1:6]..., 0.0, G.inc[1:2]...]
+	dat2las(fname, z, grd_hdr=hdr, scaleX=scaleX, scaleY=scaleY, scaleZ=scaleZ, offX=offX, offY=offY, offZ=offZ)
+end
 
-	parse_inputs_dat2las(xyz, hdr_vec)
-	n_rows, n_cols = size(xyz)
+function dat2las(fname::AbstractString, xyz; grd_hdr=Float64[], scaleX=nothing, scaleY=nothing, scaleZ=nothing, offX=nothing, offY=nothing, offZ=nothing)
 
-	if (!isempty(hdr_vec))
-		min_x, max_x, min_y, max_y, min_z, max_z = hdr_vec[1:6]
+	n_rows, n_cols = parse_inputs_dat2las(xyz, grd_hdr)
+
+	if (!isempty(grd_hdr))
+		min_x, max_x, min_y, max_y, min_z, max_z = grd_hdr[1:6]
 	end
 
 	#  Create the writer
 	laszip_writer = convert(Ptr{Ptr{Cvoid}},pointer([pointer([0])]))
-	if (laszip_create(laszip_writer) != 0)
-		msgerror(laszip_writer, "creating laszip writer")
-	end
+	(laszip_create(laszip_writer) != 0) && msgerror(laszip_writer, "creating laszip writer")
 
 	header = pointer([pointer([laszip_header()])])    # Get an empty header directly from C
 	laszip_writer = unsafe_load(laszip_writer)
@@ -43,7 +52,7 @@ function dat2las(fname::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], s
 	hdr.x_scale_factor = 1.0
 	hdr.y_scale_factor = 1.0
 	hdr.z_scale_factor = 1.0
-	if (n_cols == 3 && isempty(hdr_vec))		# The 'regular' situation
+	if (n_cols == 3 && isempty(grd_hdr))		# The 'regular' situation
 		min_x, min_y, min_z = minimum(xyz, dims=1)
 		max_x, max_y, max_z = maximum(xyz, dims=1)
 	end
@@ -56,50 +65,50 @@ function dat2las(fname::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], s
 
 	# ----------------- Find reasonable scale_factor and offset -----------------------------------------
 	if (hdr.min_x >= -360 && hdr.max_x <= 360 && hdr.min_y >= -90 && hdr.max_y <= 90)	# Assume geogs
-		hdr.x_scale_factor = isempty(scaleX) ? 1e-7 : scaleX
-		hdr.y_scale_factor = isempty(scaleY) ? 1e-7 : scaleY
+		hdr.x_scale_factor = (scaleX === nothing) ? 1e-7 : scaleX
+		hdr.y_scale_factor = (scaleY === nothing) ? 1e-7 : scaleY
 	else
-		hdr.x_scale_factor = isempty(scaleX) ? 1e-3 : scaleX
-		hdr.y_scale_factor = isempty(scaleY) ? 1e-3 : scaleY
+		hdr.x_scale_factor = (scaleX === nothing) ? 1e-3 : scaleX
+		hdr.y_scale_factor = (scaleY === nothing) ? 1e-3 : scaleY
 	end
-	hdr.z_scale_factor = isempty(scaleZ) ? 1e-2 : scaleZ
+	hdr.z_scale_factor = (scaleZ === nothing) ? 1e-2 : scaleZ
 
 	if (!isnan(hdr.min_x) && !isnan(hdr.max_x))
-		if (isempty(offX))
-			hdr.x_offset = (floor((hdr.min_x + hdr.max_x)/hdr.x_scale_factor/20000000)) * 10000000 * hdr.x_scale_factor
+		if (offX === nothing)
+			hdr.x_offset = (floor((hdr.min_x + hdr.max_x) / hdr.x_scale_factor / 20000000)) * 10000000 * hdr.x_scale_factor
 		else
 			hdr.x_offset = offX
 		end
 	end
 	if (!isnan(hdr.min_y) && !isnan(hdr.max_y))
-		if (isempty(offY))
-			hdr.y_offset = (floor((hdr.min_y + hdr.max_y)/hdr.y_scale_factor/20000000)) * 10000000 * hdr.y_scale_factor
+		if (offY === nothing)
+			hdr.y_offset = (floor((hdr.min_y + hdr.max_y) / hdr.y_scale_factor / 20000000)) * 10000000 * hdr.y_scale_factor
 		else
 			hdr.y_offset = offY
 		end
 	end
 	if (!isnan(hdr.min_z) && !isnan(hdr.max_z))
-		if (isempty(offZ))
-			hdr.z_offset = (floor((hdr.min_z + hdr.max_z)/hdr.z_scale_factor/20000000)) * 10000000 * hdr.z_scale_factor
+		if (offZ === nothing)
+			hdr.z_offset = (floor((hdr.min_z + hdr.max_z) / hdr.z_scale_factor / 20000000)) * 10000000 * hdr.z_scale_factor
 		else
 			hdr.z_offset = offZ
 		end
 	end
 	# ---------------------------------------------------------------------------------------------------
 
 	# This is the case where we are storing a grid pretending it's a regular LAZ file. Must hijack some header members
-	if (!isempty(hdr_vec))
+	if (!isempty(grd_hdr))
 		hdr.x_scale_factor = hdr.z_scale_factor		# Because in fact we only have zz's
 		hdr.y_scale_factor = hdr.z_scale_factor
 		hdr.x_offset = hdr.z_offset
 		hdr.y_offset = hdr.z_offset
 		hdr.number_of_point_records = UInt32(ceil(n_rows / 3))
 		hdr.global_encoding = 32768;	# Use this number to codify as GRID. bin(UInt16(32768)) = "1000000000000000"
 
-		one = (hdr_vec[7] == 0 ? 1 : 0)
-		hdr.project_ID_GUID_data_1 = hdr_vec[7]
-		hdr.project_ID_GUID_data_2 = round(UInt16, (hdr.max_y - hdr.min_y) / hdr_vec[8]) + one 	# n_rows in 2D array
-		hdr.project_ID_GUID_data_3 = round(UInt16, (hdr.max_x - hdr.min_x) / hdr_vec[9]) + one 	# n_cols in 2D array
+		one = (grd_hdr[7] == 0 ? 1 : 0)
+		hdr.project_ID_GUID_data_1 = UInt32(grd_hdr[7])
+		hdr.project_ID_GUID_data_2 = round(UInt16, (hdr.max_y - hdr.min_y) / grd_hdr[8]) + one 	# n_rows in 2D array
+		hdr.project_ID_GUID_data_3 = round(UInt16, (hdr.max_x - hdr.min_x) / grd_hdr[9]) + one 	# n_cols in 2D array
 	end
 
 	# Save back the header to its C pointer
@@ -167,25 +176,22 @@ function dat2las(fname::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], s
 	end
 
 	# Close the writer
-	if (laszip_close_writer(laszip_writer) != 0)
-		msgerror(laszip_writer, "closing laszip writer")
-	end
+	(laszip_close_writer(laszip_writer) != 0) && msgerror(laszip_writer, "closing laszip writer")
+
 	# Destroy the writer
-	if (laszip_destroy(laszip_writer) != 0)
-		msgerror(laszip_writer, "destroying laszip writer")
-	end
+	(laszip_destroy(laszip_writer) != 0) && msgerror(laszip_writer, "destroying laszip writer")
 
+	return nothing
 end
 
 # --------------------------------------------------------------------------
-function parse_inputs_dat2las(xyz, hdr_vec)
+function parse_inputs_dat2las(xyz, grd_hdr)
 # Check validity of input and in future will parse string options
 
-	n_rows, n_cols = size(xyz)
+	(!isempty(grd_hdr) && length(grd_hdr) < 7) && error("HDR argument does not have at least 7 elements")
+	n_rows, n_cols = isa(xyz, Matrix) ? size(xyz) : (length(xyz), 1)
 	((n_cols != 3 && n_cols != 1)) && error("Input array can only have 1 or 3 columns OR be a 2D array")
-	if (!isempty(hdr_vec))
-		(length(hdr_vec) < 7) && error("HDR argument does not have at least 7 elements")
-	end
+	return n_rows, n_cols
 end
 
 const xyz2laz  = dat2las			# Alias