diff --git a/src/common_options.jl b/src/common_options.jl index 10c24a101..fade5ddc3 100644 --- a/src/common_options.jl +++ b/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 diff --git a/src/gmt_main.jl b/src/gmt_main.jl index 90b78156a..930d36694 100644 --- a/src/gmt_main.jl +++ b/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 diff --git a/src/imgtiles.jl b/src/imgtiles.jl index c60bde956..fdf8e1307 100644 --- a/src/imgtiles.jl +++ b/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 diff --git a/src/webmapserver.jl b/src/webmapserver.jl index 300ec88fb..877f625d6 100644 --- a/src/webmapserver.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 3dacfcfa0..ceec0eb9f 100644 --- a/test/runtests.jl +++ b/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