Let the wmsreader fully deal with limits given by tile names.

src/common_options.jl

@@ -176,12 +176,13 @@ function parse_R(d::Dict, cmd::String, O::Bool=false, del::Bool=true, RIr::Bool=
 			CTRL.limits[7:7+length(limits)-1] = limits		# The plot limits
 			(opt_R != " -Rtight" && opt_R !== nothing && limits != zeros(4) && all(CTRL.limits[1:4] .== 0)) &&
 				(CTRL.limits[1:length(limits)] = limits)	# And this makes data = plot limits, IF data is empty.
-			if (count_chars(opt_R, ',') == 2)				# A XYZ tile address (like "7829,6374,14")
+			if (istilename(opt_R))							# A XYZ or quadtree tile address (like "7829,6374,14")
 				opt_R = sprintf(" -R%.12g/%.12g/%.12g/%.12g", limits[1], limits[2], limits[3], limits[4])
 			elseif contains(opt_R, "+r")
 				CTRL.limits[13] = 1.0						# To know that -R...+r was used.
 			end
-		catch
+		catch err
+			@warn err
 			CTRL.limits .= 0.0
 			IamModern[1] = bak
 		end
@@ -363,8 +364,12 @@ function opt_R2num(opt_R::String)::Vector{Float64}
 		for k = 2:lastindex(rs)  limits[k] = parse(Float64, rs[k])  end
 		#if (isdiag)  limits[2], limits[4] = limits[4], limits[2]  end
 		# Don't know anymore how -R...+r limits should be stored in CTRL.limits
-	elseif (count_chars(opt_R, ',') == 2)		# A XYZ tile address
-		limits = mosaic(scan_opt(opt_R, "-R"), mesh=true).bbox
+	elseif (istilename(opt_R))						# A XYZ or quadtree tile address
+		t = scan_opt(opt_R, "-R")
+		limits = mosaic(t, mesh=true).bbox
+		zl::Int = contains(t, ",") ? parse(Int, t[findlast(',', t)+1:end]) : length(t)
+		inc = (360 / 2 ^ zl) / 256					# Increment in degrees at this zoom level
+		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.
@@ -4159,7 +4164,7 @@ function showfig(d::Dict, fname_ps::String, fname_ext::String, opt_T::String, K:
 	end
 	CTRL.limits .= 0.0;		CTRL.figsize .= 0.0;	CTRL.proj_linear[1] = true;		# Reset these for safety
 	CTRL.pocket_J[1], CTRL.pocket_J[2], CTRL.pocket_J[3], CTRL.pocket_J[4] = "", "", "", "   ";
-	CTRL.pocket_R[1] = "";	isJupyter[1] = false;	CTRL.pocket_call .= nothing
+	CTRL.pocket_R[1:2] .= "";	isJupyter[1] = false;	CTRL.pocket_call .= nothing
 	CURRENT_VIEW[1] = ""
 	return retPluto ? WrapperPluto(out) : nothing	# retPluto should make it all way down to base so that Plut displays it
 end

src/gmt_main.jl

@@ -1401,7 +1401,7 @@ function resetGMT(dorestart::Bool=true)
 	DEF_FIG_AXES[1] = DEF_FIG_AXES_BAK;		DEF_FIG_AXES3[1] = DEF_FIG_AXES3_BAK;
 	CTRL.pocket_J[1], CTRL.pocket_J[2], CTRL.pocket_J[3], CTRL.pocket_J[4] = "", "", "", "   ";
 	CTRL.IamInPaperMode[:] = [false, true];	IamInset[1], IamInset[2] = false, false
-	CTRL.pocket_call[1:3] .= nothing;	CTRL.pocket_R[1] = "";	CTRL.figsize .= 0.0
+	CTRL.pocket_call[1:3] .= nothing;	CTRL.pocket_R[1:2] .= "";	CTRL.figsize .= 0.0
 	CTRL.XYlabels[1] = "";	CTRL.XYlabels[2] = "";	CTRL.returnPS[1] = false
 	(dorestart) && (gmt_restart(); clear_sessions())
 end

src/imgtiles.jl

@@ -163,6 +163,7 @@ viz(I, coast=true)
 """
 function mosaic(; pt_radius=6378137.0, provider="", zoom::Int=0, cache::String="",
                 mapwidth=15, dpi=96, date::String="", verbose::Int=0, kw...)
+	isempty(kw) && return mosaic(zoom)		# Call the method that only prints the zoom levels table.
 	d = KW(kw)
 	((opt_R = parse_R(d, "")[1]) == "") && error("To use the 'mosaic' function without the 'lon & lat' arguments you need to specify the 'region' option.")
 	ll = opt_R2num(opt_R)
@@ -171,7 +172,7 @@ function mosaic(; pt_radius=6378137.0, provider="", zoom::Int=0, cache::String="
            dpi=dpi, date=date, verbose=verbose, d...)
 end
 
-# This methos is mostly for calls from python's juliacall that used PyList (because dumb Py consider this a list: [1.0, 2.6])
+# This method is mostly for calls from python's juliacall that used PyList (because dumb Py consider this a list: [1.0, 2.6])
 function mosaic(lon::AbstractVecOrMat, lat::AbstractVecOrMat; pt_radius=6378137.0, provider="", zoom::Int=0, cache::String="",
                 mapwidth=15, dpi=96, verbose::Int=0, date::String="", key::String="", kw...)
 	_lon::Vector{Float64}, _lat::Vector{Float64} = vec(Float64.(lon)), vec(Float64.(lat))
@@ -443,6 +444,55 @@ function mosaic(lon::Vector{<:Float64}, lat::Vector{<:Float64}; pt_radius=637813
 	return I
 end
 
+# ---------------------------------------------------------------------------------------------------
+"""
+    mosaic([zoom::Int=??])
+
+Print a table with the zoom level characteristics in terms of tile sizes, resolutions, typical use.
+
+If the `zoom` option is used then the table is printed with that zoom level only.
+
+# Example
+```jldoctest
+julia> mosaic(zoom=10)
+┌───────┬────────────────┬────────────┬────────────────┬────────────────────┐
+│ Level │     Tile width │  m / pixel │         ~Scale │        Examples of │
+│       │ ° of longitude │ on Equator │                │ areas to represent │
+├───────┼────────────────┼────────────┼────────────────┼────────────────────┤
+│    10 │          0.352 │    153.054 │ 1:500 thousand │  metropolitan area │
+└───────┴────────────────┴────────────┴────────────────┴────────────────────┘
+```
+"""
+function mosaic(zoom)
+	(zoom < 1 || zoom > 22) && (zoom = 0)
+	hdr = (["Level", "Tile width", "m / pixel", "~Scale", "Examples of"], ["", "° of longitude", "on Equator", "", "areas to represent"])
+	data = Any[
+			1  180 78_272 "1:250 million" ""
+			2   90 39_136 "1:150 million" "subcontinental area"
+			3   45 19_568 "1:70 million"  "largest country"
+			4   22.5  9_784 "1:35 million" ""
+			5 	11.25 4_892 "1:15 million" "large African country"
+			6 	5.625 2_446 "1:10 million" "large European country"
+			7 	2.813 1_223 "1:4 million" "small country, US state"
+			8 	1.406 611.388 "1:2 million" ""
+			9 	0.703 305.694 "1:1 million" "wide area, large metropolitan area"
+			10 	0.352 153.054 "1:500 thousand" "metropolitan area"
+			11 	0.176 76.532 "1:250 thousand" "city"
+			12 	0.088 38.266 "1:150 thousand" "town, or city district"
+			13 	0.044 19.133 "1:70 thousand" "village, or suburb"
+			14 	0.022 9.566 "1:35 thousand" ""
+			15 	0.011 4.783 "1:15 thousand" "small road"
+			16 	0.005 2.174 "1:8 thousand" "street"
+			17 	0.003 1.305 "1:4 thousand" "block, park, addresses"
+			18 	0.001 0.435 "1:2 thousand" "some buildings, trees"
+			19 	0.0005 0.217 "1:1 thousand" "local highway and crossing details"
+			20 	0.00025 0.109 "1:5 hundred" "A mid-sized building"
+			21 	0.000125 0.054 "1:250" ""
+		]
+	println("\t\t\t\tZoom levels")
+	(zoom == 0) ? pretty_table(data; header=hdr) : pretty_table(data[zoom:zoom, :]; header=hdr)
+end
+
 # ---------------------------------------------------------------------------------------------------
 """
     getprovider(name, zoom::Int; variant="", date::String="", key::String="")
@@ -924,3 +974,17 @@ function geocoder(address::String; options=String[])
 	D.ds_bbox = [BB[3], BB[4], BB[1], BB[2]]
 	return D
 end
+
+# ------------------------------------------------------------------------------------------------
+"""
+    istilename(s::AbstractString)
+
+Check if the string `s` is a XYZ or quadtree tile name. Useful for parse_R() and others that can
+than extract the tile limits and associated resolution.
+"""
+function istilename(s::AbstractString)::Bool
+	(count_chars(s, ',') == 2) && return true
+	contains(s, "-R") && occursin(r"^[0-3]+$", s[findfirst('R', s)+1:end]) && return true	# Accepts also that 's' is an opt_R
+	occursin(r"^[0-3]+$", s) && return true
+	return false
+end

src/webmapserver.jl

@@ -105,12 +105,26 @@ Read the `layer` number provided by the service from which the `wms` type was cr
    both of these forms are allowed: `layer=3` or `layer="Invented layer name"`
 
 ### `kwargs` is the keywords/values pairs used to set
-- `region | limits`: The region limits. This can be a Tuple or Array with 4 elements defining the `(xmin, xmax, ymin, ymax)`
-   or a string defining the limits in all ways that GMT can recognize. When the layer has the data projected, we can
-   a Tuple or Array with 3 elements `(lon0, lat0, width)`, where first two set the center of a square in geographical
-   coordinates and the third (`width`) is the width of that box in meters.
+- `region | limits`: The region limits. Options are:
+   - A Tuple or Array with 4 elements defining the `(xmin, xmax, ymin, ymax)` or a string defining the
+     limits in all ways that GMT can recognize. 
+   - When the layer has the data projected, we can pass a Tuple or Array with 3 elements `(lon0, lat0, width)`,
+     where first two set the center of a square in geographical coordinates and the third (`width`) is the
+     width of that box in meters.
+   - A ``mosaic`` tile name in the form ``X,Y,Z`` or a quadtree. Example: ``region="7829,6374,14"``. See the ``mosaic``
+     function manual for more information. This form also sets the default cellsize for that tile. NOTE:
+     this is a geographical coordinates only that implicitly sets ``geog=true``. See below on how to change
+     the default resolution.
 - `cellsize | pixelsize | resolution | res`: Sets the requested cell size in meters [default]. Use a string appended with a 'd'
    (e.g. `resolution="0.001d"`) if the resolution is given in degrees. This way works only when the layer is in geogs.
+- `zoom or refine`: When the region is derived from a ``mosaic`` tile name, the default is to get an image with 256 columns
+   and _n_ rows where _n_ depends on latitude. So, either the area is large and consequently the resolution is low, or
+   the inverse (small area and resolution is high). To change this status, use the `zoom` or `refine` options.
+   - `zoom`: an integer >= 1 that for each increment duplicates the base resolution by 2. _e.g._, `zoom=2`
+      quadruplicates the default resolution. This option is almost redundant with the `refine`, but is offered
+      for consistency with the ``mosaic`` function.
+   - `refine`: an integer >= 2 multiplication factor that is used to increment the resolution by factor. _e.g._, `refine=2`
+      duplicates the image resolution.
 - `size`: Alternatively to the the `cellsize` use this option, a tuple or array with two elements, to specify
    the image dimensions. Example, `size=(1200, 100)` to get an image with 1200 rows and 100 columns.
 - `time`: Some services provide data along time. Use this option to provide a time string as provided by DateTime.
@@ -132,7 +146,7 @@ A GMTimage
     img = wmsread(wms, layer="MODIS_Terra_CorrectedReflectance_TrueColor", region=(9,22,32,43), time="2021-10-29T00:00:00", pixelsize=750);
     imshow(img, proj=:guess)
 """
-function wmsread(wms::WMS; layer=0, time::String="", kw...)
+function wmsread(wms::WMS; layer=0, time::String="", kw...)::GMTimage
 	layer_n = get_layer_number(wms, layer)
 	str, dim_x, dim_y = wms_helper(wms; layer=layer_n, time=time, kw...)
 	opts = ["-outsize", "$(dim_x)", "$(dim_y)"]
@@ -212,12 +226,11 @@ function wms_helper(wms::WMS; layer=0, kw...)
 			half_w::Float64 = reg[3] / 2
 			lims = [center[1]-half_w, center[1]+half_w, center[2]-half_w, center[2]+half_w]
 		end
-		BB = @sprintf("%.12g,%.12g,%.12g,%.12g", lims[1], lims[3], lims[2], lims[4])
 	else
-		_, opt_R = parse_R(d, "")
-		lims = opt_R2num(opt_R)
-		BB = @sprintf("%.12g,%.12g,%.12g,%.12g", lims[1], lims[3], lims[2], lims[4])
+		parse_R(d, "")			# This fills CTRL.limits, which is what we need here
+		lims = CTRL.limits[1:4]
 	end
+	BB = @sprintf("%.12g,%.12g,%.12g,%.12g", lims[1], lims[3], lims[2], lims[4])
 	(wms.layer[layer_n].bbox[1] > lims[1] || wms.layer[layer_n].bbox[2] < lims[2] || wms.layer[layer_n].bbox[3] > lims[3] || wms.layer[layer_n].bbox[4] < lims[4]) && @warn("Requested region overflows this layer BoundingBox")
 
 	function loc_getsize(res::Float64, lims::Vector{<:Real})::Tuple{Int, Int}
@@ -245,6 +258,19 @@ function wms_helper(wms::WMS; layer=0, kw...)
 			end
 			dim_x, dim_y = loc_getsize(_res, lims)
 		end
+	elseif (CTRL.limits[2] != "")		# Indirect check if -R was set via a mosaic tile name. A bit FRAGILE
+		try _res = tryparse(Float64, CTRL.pocket_R[2])
+		catch; error("Programming error. 'CTRL.pocket_R[2]' is empty")
+		end
+		if ((fact_ = find_in_dict(d, [:refine])[1]) !== nothing)
+			((fact = tryparse(Int, fact_)) === nothing) && error("The 'refine' option must be an integer")
+			_res /= fact
+		elseif ((zoom = find_in_dict(d, [:zoom])[1]) !== nothing)
+			@assert (isa(zoom, Int) && zoom > 0) "'zoom' is not an integer > 0"
+			_res /= 2 ^ zoom
+		end
+		dim_x, dim_y = loc_getsize(_res, lims)
+		SRS = "EPSG:4326"			# This is a geogs case only
 	else
 		error("Must provide either the 'size' or the 'resolution' option")
 	end

test/runtests.jl

@@ -38,6 +38,7 @@ using Dates, Printf#, Logging
 		GMT.wmstest(wms, layer=33, region=(iso="PT"), res=100);
 		GMT.wmstest(wms, layer=37, region=(-8,-7,38,39), res="0.001d")
 		GMT.wmstest(wms, layer=37, region=(-8,-7,38,39), res=100)
+		@test GMT.wmstest(wms, layer=37, region="7829,6374,14", zoom=3, size=true) == (1635, 2048)
 	catch
 	end