Call grdinfo -C, when possible, via the gmt_grdinfo_C() function. (#1869)

Author: joa-quim

src/GMT.jl

@@ -163,6 +163,7 @@ export
 	ogr2ogr, gdalrasterize, gdalbuildvrt, readgeom, readraster, setgeotransform!, setnodata!, setproj!, destroy,
 	delaunay, dither, buffer, centroid, intersection, intersects, polyunion, overlaps, fromWKT, fillnodata!, fillnodata,
 	concavehull, convexhull, difference, symdifference, distance, geodesicarea, geomarea, pointalongline, polygonize, simplify,
+	boundary, crosses, disjoint, envelope, envelope3d, equals, touches, uniqueind, within, 
 
 	wkbUnknown, wkbPoint, wkbPointZ, wkbLineString, wkbLineStringZ, wkbPolygon, wkbPolygonZM, wkbMultiPoint, wkbMultiPointZ,
 	wkbMultiLineString, wkbMultiPolygon, wkbGeometryCollection, wkbPoint25D, wkbLineString25D, wkbPolygon25D, wkbMultiPoint25D,
@@ -197,8 +198,7 @@ export
 	gunique, sortslicesperm,
 	hampel, hampel!, lowess, whittaker,
 
-	Ginnerjoin, Gouterjoin, Gleftjoin, Grightjoin, Gcrossjoin, Gsemijoin, Gantijoin, spatialjoin,
-	groupby, stats,
+	Ginnerjoin, Gouterjoin, Gleftjoin, Grightjoin, Gcrossjoin, Gsemijoin, Gantijoin, spatialjoin, groupby, stats,
 
 	anaglyph,
 	lazinfo, lazread, lazwrite, lasread, laswrite,
@@ -424,6 +424,7 @@ using .Laszip
 	grdview(rand(Float32,32,32), Vd=2);
 	grdinfo(mat2grid(rand(Float32,4,4)));
 	Glix=gmt("grdmath", "-R0/10/0/10 -I2 X");
+	gmt_grdinfo_C(Glix);
 	grdcontour(Glix);
 	grd2cpt(Glix);
 	grd2xyz(Glix);

src/common_options.jl

@@ -176,7 +176,7 @@ function parse_R(d::Dict, cmd::String; O::Bool=false, del::Bool=true, RIr::Bool=
 			opt_r = parse_r(d, "")[2]
 			(opt_r == "") && (cmd *= " -r" * ((val.registration == 0) ? "g" : "p"))
 		elseif (!isempty(val_str) && ((t = guess_T_from_ext(val_str)) == " -Tg" || t == " -Ti"))		# A file name
-			info::Matrix{Float64} = gmt("grdinfo -C " * val_str).data			# Get the grid info
+			info::Matrix{Float64} = gmt_grdinfo_C(val_str).data			# Get the grid info
 			opt_R = @sprintf(" -R%.12g/%.12g/%.12g/%.12g", info[1,1], info[1,2], info[1,3], info[1,4])
 			(opt_I = parse_I(d, "", [:I :inc :increment :spacing], "I", true)) == "" &&
 				(opt_I = @sprintf(" -I%.8g/%.8g", info[7], info[8]))	# The grid increment
@@ -388,8 +388,8 @@ function opt_R2num(opt_R::String)::Vector{Float64}
 		CTRL.pocket_R[2] = "$inc"
 	elseif (opt_R != " -R" && opt_R != " -Rtight")	# One of those complicated -R forms. Ask GMT the limits (but slow. It takes 0.2 s)
 		# If opt_R is not a grid's name, we are f.
-		((ind = findfirst("-R@", opt_R)) !== nothing) && return gmt("grdinfo " * opt_R[ind[3]:end] * " -C")[1:4]	# should be a cache file
-		(((f = guess_T_from_ext(opt_R)) == " -Tg") || f == " -Ti") && return gmt("grdinfo " * opt_R * " -C")[1:4]	# any local file
+		((ind = findfirst("-R@", opt_R)) !== nothing) && return gmt_grdinfo_C(opt_R[ind[3]:end])[1:4]	# should be a cache file
+		(((f = guess_T_from_ext(opt_R)) == " -Tg") || f == " -Ti") && return gmt_grdinfo_C(opt_R)[1:4]	# any local file
 
 		(opt_R == " -R=WD") && return [-180.0, 180., -90., 90.]		# World map. Shit is if [0 360] is wanted.
 		kml::GMTdataset{Float64, 2} = gmt("gmt2kml " * opt_R, [0 0])	# for example, opt_R = " -RPT"
@@ -2740,7 +2740,7 @@ function get_cpt_set_R(d::Dict, cmd0::String, cmd::String, opt_R::String, got_fn
 	if (isa(arg1, GItype))			# GMT bug, -R will not be stored in gmt.history
 		range::Vector{Float64} = isa(arg1.range, Float64) ? vec(arg1.range) : vec(arg1.range[1:6])	# Beyond 6 can be a Time
 	elseif (cmd0 != "" && cmd0[1] != '@')
-		info = grdinfo(cmd0 * " -C");	range = vec(info.data)
+		info = gmt_grdinfo_C(cmd0);	range = vec(info.data)
 	end
 
 	cpt_opt_T = ""
@@ -4637,8 +4637,7 @@ function regiongeog(GI::GItype)::Tuple
 end
 function regiongeog(fname::String)::Tuple
 	((prj = getproj(fname, wkt=true)) == "") && (@warn("Input grid/image has no projection info"); return ())
-	#info = grdinfo(fname, C=true);		# It should also report the
-	info = gmt("grdinfo -C " * fname);
+	info = gmt_grdinfo_C(fname)
 	c = xy2lonlat([info.data[1] info.data[3]; info.data[2] info.data[4]]; s_srs=prj)
 	tuple(c...)
 end

src/contourf.jl

@@ -151,7 +151,7 @@ function contourf(cmd0::String, arg1, arg2, first::Bool, d::Dict{Symbol,Any})
 		#isa(arg2, GMTcpt) ? d[:N] = arg2 : (isa(arg1, GMTcpt) ? d[:N] = arg1 : d[:N] = true)
 		if (isempty(CPT) && isempty(CPT_arg))
 			if (cmd0 != "")
-				info = grdinfo(cmd0 * " -C")
+				info = gmt_grdinfo_C(cmd0)
 				C_inc, min, max = gen_contour_vals(info.data[5:6], C_int)
 			else
 				C_inc, min, max = gen_contour_vals(arg1, C_int)

src/gdal/gdal_tools.jl

@@ -253,7 +253,7 @@ function gdaldem(indata, method::String, opts::Vector{String}=String[]; dest="/v
 		append!(opts, ["-compute_edges", "-b", band])
 		if ((val = find_in_dict(d, [:scale])[1]) === nothing)
 			if (isa(indata, GMTgrid) && (occursin("longlat", indata.proj4) || occursin("latlong", indata.proj4)) ||
-			                             gmt("grdinfo -Cn", indata).data[end] == 1)
+			                             GMT.gmt_grdinfo_C(indata).data[end] == 1)
 				append!(opts, ["-s", "111120"])
 			end
 		else

src/gdal_utils.jl

@@ -488,7 +488,7 @@ by `grdinfo`. An unfortunate consequence of using either GMT or GDAL to inquire
 `x_inc, y_inc` when reading with GDAL, and `z_min, z_max` when reading with GMT.
 """
 function getregion(input; pad=0, xSize=0, ySize=0, gridreg::Bool=false, sds::Int=0, GMT=false)::NTuple{6, Float64}
-	(GMT == 1 && isa(input, String)) && return vec(grdinfo(input, C=:n).data)
+	(GMT == 1 && isa(input, String)) && return vec(gmt_grdinfo_C(input).data)
 	dataset = gd2gmt_helper(input, sds; short=true)[1]
 	(xSize == 0) && (xSize = Gdal.width(dataset))
 	(ySize == 0) && (ySize = Gdal.height(dataset))

src/gmt_main.jl

@@ -236,6 +236,8 @@ end
 gmt_GMTcpt(cmd::String, args...)::Union{GMTcpt, Nothing} = gmt(cmd, args...)::Union{GMTcpt, Nothing}
 gmt_GMTdataset(cmd::String, args...)::Union{GMTdataset{Float64,2}, Vector{<:GMTdataset{Float64,2}}} = gmt(cmd, args...)
 gmt_GMTgrid(cmd::String, args...)::GMTgrid{Float32,2} = gmt(cmd, args...)::GMTgrid{Float32,2}
+gmt_grdinfo_C(name::String)::GMTdataset{Float64, 2} = gmt("grdinfo -C " * name)::GMTdataset{Float64, 2}
+gmt_grdinfo_C(GI::GItype)::GMTdataset{Float64, 2} = gmt("grdinfo -C ", GI)::GMTdataset{Float64, 2}
 
 # -----------------------------------------------------------------------------------------------
 function gmt_restart(restart::Bool=true)

src/grdinfo.jl

@@ -80,6 +80,3 @@ function grdinfo_helper(cmd0::String, arg1; kwargs...)::Union{GMTdataset, String
 	end
 	return R
 end
-
-# ---------------------------------------------------------------------------------------------------
-#grdinfo(arg1; kw...) = grdinfo("", arg1; kw...)

src/grdinterpolate.jl

@@ -122,7 +122,7 @@ function grdinterp_helper(cmd0::String, arg1; allcols::Bool=false, gdal=false, k
 		common_grd(d, "grdinterpolate " * cmd, arg1..., arg2)
 	else
 		if (cmd0 != "" && !contains(cmd, " -R") && val !== nothing)		# Since GMT will read layer-by-layer, better to limit to the region framing the points
-			D = grdinfo(cmd0, C=true)
+			D = gmt_grdinfo_C(cmd0)
 			x_min = max(D[1], pts.ds_bbox[1]-D[9]*3);	x_max = min(D[2], pts.ds_bbox[2]+D[9]*3)
 			y_min = max(D[3], pts.ds_bbox[3]-D[10]*3);	y_max = min(D[4], pts.ds_bbox[4]+D[10]*3)
 			cmd *= " -R$x_min/$x_max/$y_min/$y_max"

src/grdvector.jl

@@ -178,8 +178,8 @@ function parse_Q_grdvec(d::Dict, symbs::Array{<:Symbol}, len::String="", stop::S
 end
 
 # ---------------------------------------------------------------------------------------------------
-get_grdinfo(grd::String,  opt_R::String) = gmt("grdinfo -C" * opt_R * " " * grd).data
-get_grdinfo(grd::GMTgrid, opt_R::String) = gmt("grdinfo -C" * opt_R, grd).data
+get_grdinfo(grd::String,  opt_R::String) = gmt_grdinfo_C(opt_R * " " * grd).data
+get_grdinfo(grd::GMTgrid, opt_R::String)::Matrix{Float64} = gmt("grdinfo -C" * opt_R, grd).data
 
 # ---------------------------------------------------------------------------------------------------
 grdvector!(arg1, arg2; kw...) = grdvector(arg1, arg2; first=false, kw...)

src/imshow.jl

@@ -251,7 +251,7 @@ function snif_GI_set_CTRLlimits(G_I)::Bool
 	(isa(G_I, String) && (G_I[1] == '@' || startswith(G_I, "http"))) && return false	# Remotes are very dangerous to sniff in
 
 	# Do not call grdinfo over grid/images already in memory
-	range::Vector{Float64} = isa(G_I, String) ? vec(grdinfo(G_I, C=:n).data) : G_I.range[1:6]
+	range::Vector{Float64} = isa(G_I, String) ? vec(gmt_grdinfo_C(G_I).data) : G_I.range[1:6]
 	if ((isa(G_I, String) && range[2] != range[9] && range[4] != range[10]) ||
 		(!isa(G_I, String) && range[2] != size(G_I, 2) && range[4] != size(G_I, 1)))
 		CTRL.limits[1:4], CTRL.limits[7:10] = range[1:4], range[1:4]