Merge pull request #2000 from GenericMappingTools/mosaic-ortosat23

Let mosaic() also deal with the Ortosat23 Portuguese imagerie via new provider.

Author: joa-quim

src/extras/webmapserver.jl

@@ -41,7 +41,7 @@ returned by the first form.
 """
 function wmsinfo(server::String)::WMS
 	w = gdalinfo("WMS:" * server)
-	(w === nothing) && error("unable to open '$(server)'")
+	(w === nothing || w === "") && error("unable to open '$(server)'")
 	sds = split(w, "SUBDATASET_")
 	lastguy = collect(2:2:length(sds))[end]
 

src/gdal_utils.jl

@@ -1067,7 +1067,10 @@ function lonlat2xy(xy::Matrix{T}, t_srs_=nothing; t_srs=nothing, s_srs=prj4WGS84
 	opts = (size(xy,2) == 2) ? ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite"] :
 	                           ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite", "-dim", "XYZ"]
 	D = ogr2ogr(xy, opts)
-	return D.data		# Return only the array because that's what was sent in
+	# ogr2ogr/gd2gmt returns a single GMTdataset when the input has one point, but a
+	# Vector{GMTdataset} when it has several (each row becomes a separate point feature).
+	# Aggregate back into a single matrix matching the input shape.
+	return isa(D, Vector) ? reduce(vcat, (d.data for d in D)) : D.data
 end
 
 lonlat2xy(D::GMTdataset, t_srs_=nothing; t_srs=nothing, s_srs=prj4WGS84) = lonlat2xy([D], t_srs_; t_srs=t_srs, s_srs=s_srs)
@@ -1118,7 +1121,10 @@ function xy2lonlat(xy::Matrix{<:Real}, s_srs_=""; s_srs="", t_srs=prj4WGS84)
 	(s_srs == "") && error("Must specify at least the source referencing system.")
 	D = ogr2ogr(xy, ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite"])
 	isa(D, Vector) && contains(s_srs, "=lee_os") && length(D) > 1 && return [D[1][1,1] D[1][1,2]; D[2][2,1] D[2][2,2]]	# GDAL BUG?
-	return D.data		# Return only the array because that's what was sent in
+	# ogr2ogr/gd2gmt returns a single GMTdataset when the input has one point, but a
+	# Vector{GMTdataset} when it has several (each row becomes a separate point feature).
+	# Aggregate back into a single matrix matching the input shape.
+	return isa(D, Vector) ? reduce(vcat, (d.data for d in D)) : D.data
 end
 
 xy2lonlat(D::GMTdataset, s_srs_=""; s_srs="", t_srs=prj4WGS84) = xy2lonlat([D], s_srs_; s_srs=s_srs, t_srs=t_srs)

src/imgtiles.jl

@@ -135,6 +135,16 @@ Get image tiles from a web map tiles provider for given longitude, latitude coor
   - "Bing" (the default), "Google", "OSM", "Esri" or a custom provider.
   - A `Provider` type from the ``TileProviders.jl`` package. You must consult the documentation of that package
     for more details on how to choose a *provider*.
+  - "ortosat23" / "Ortosat23": Portuguese DGT orthophotos (2023) served via WMS. Always returned in
+    WGS84 geographical coordinates. The spelling controls the output pipeline:
+      * `provider="ortosat23"` (lowercase): single-warp **approximate** path. Cheap (one warp). Output
+        size differs slightly from the same-zoom Bing/Google result because the WMS source SRS
+        (EPSG:3763) and Bing's source SRS (EPSG:3857) have different aspect after warping to geog.
+      * `provider="Ortosat23"` (capital `O`): double-warp **exact** path. Two `gdalwarp` calls
+        (3763 -> 3857 -> latlong) replicate Bing's pipeline so the output bbox and pixel count match
+        Bing exactly at the same `zoom` / `neighbors`. More expensive (two warps over the full image).
+    For native EPSG:3763 metric coordinates or custom pixel resolution / `zoom` / `refine` factors,
+    call `wmsread` directly (see `? wmsread`).
 - `zoom`: Zoom level (0 for automatic). A number between 0 and ~19. The maximum is provider and area dependent.
   If `zoom=0`, the zoom level is computed automatically based on the `mapwidth` and `dpi` options.
 - `cache`: Full name of the the cache directory where to save the downloaded tiles. If empty, a cache
@@ -482,6 +492,12 @@ function mosaic(lon::Vector{<:Float64}, lat::Vector{<:Float64}; pt_radius=637813
 	# Return here if user wants a GMTdataset with the coordinates of the tiles
 	quadTiles && return quadbounds(quad_; quadkey=quadkey)[1]
 
+	# -------- WMS-backed providers (ortosat23): fetch one image covering the tile-matrix bbox ---------
+	if (provider_code == "ortosat23" || provider_code == "Ortosat23")
+		return _mosaic_ortosat23(provider_url, variant, lon_mm, lat_mm, mMo, nMo, zoom, pt_radius;
+		                        exact = (provider_code == "Ortosat23"))
+	end
+
 	# ------------------ Get the tiles and build up the image --------------------
 	cache, cache_supp = completeCacheName(cache, zoom, provider_code; variant=variant)
 	if !isa(tile_url, Matrix)
@@ -513,6 +529,63 @@ function mosaic(lon::Vector{<:Float64}, lat::Vector{<:Float64}; pt_radius=637813
 	return I
 end
 
+# ---------------------------------------------------------------------------------------------------
+# Helper for the 'ortosat23' / 'Ortosat23' providers: instead of fetching individual XYZ tiles, query
+# the DGT WMS service once with the bounding box covering the tile-matrix area. Returns a GMTimage in
+# geographical coordinates (WGS84). The `variant` selects the WMS layer ("CorVerdadeira" or "FalsaCor");
+# when no exact/partial match is found, layer 2 is used.
+#
+# Two output pipelines, selected by the `exact` flag:
+#  - `exact=false` (provider="ortosat23"): single-warp approximate path. Passes the geog bbox to
+#    `wmsread` with a `pixelsize` derived from the Web Mercator ground resolution at the bbox
+#    mid-latitude. `wms_helper` reprojects bbox to EPSG:3763, gdaltranslate fetches the image and the
+#    internal gdalwarp warps to latlong with auto-sizing. Cheap (one warp), but output size differs
+#    slightly from the same-zoom Bing/Google result because the source SRS (3763) and Bing's source
+#    SRS (3857) have different aspect after warping to geog.
+#  - `exact=true` (provider="Ortosat23"): double-warp exact path. Projects the bbox to EPSG:3763
+#    ourselves and passes those *projected* lims plus `size=(mMo*256, nMo*256)` to `wmsread`. Since
+#    the lims are no longer in geog range, `wms_helper` skips its internal warp and the image stays
+#    in 3763. We then warp 3763 -> 3857 (with `-ts` and `-te` matching Bing's intermediate mosaic in
+#    spherical Mercator) and finally 3857 -> latlong (auto-size, identical to Bing's last step). The
+#    output bbox/size then matches Bing exactly. Two warps over the full image.
+function _mosaic_ortosat23(url::String, variant::String, lon_mm::Vector{Float64},
+                           lat_mm::Vector{Float64}, mMo::Int, nMo::Int, zoom::Int,
+                           pt_radius::Float64; exact::Bool=false)::GMTimage
+	wms = wmsinfo(url)
+	layer_n = 2
+	if (variant != "")
+		v = lowercase(variant)
+		idx = findfirst(n -> occursin(v, lowercase(n)), wms.layernames)
+		(idx !== nothing) && (layer_n = idx)
+	end
+
+	if (!exact)				# Single-warp approximate path
+		bbox = (minimum(lon_mm), maximum(lon_mm), minimum(lat_mm), maximum(lat_mm))
+		lat_mid = (bbox[3] + bbox[4]) / 2
+		res_m = 156543.034 * cosd(lat_mid) / 2^zoom	# Web Mercator ground resolution at zoom & lat_mid
+		return wmsread(wms; layer=layer_n, region=bbox, pixelsize=res_m)
+	end
+
+	# --- Double-warp exact path ---
+	# 1) Project geog bbox -> layer SRS (EPSG:3763) ourselves so wmsread keeps the image in 3763.
+	srs_str::String = wms.layer[layer_n].srs
+	t_srs_3763::String = epsg2proj(parse(Int, srs_str[6:end]))
+	corners_3763::Matrix{Float64} = lonlat2xy([lon_mm[1] lat_mm[1]; lon_mm[2] lat_mm[2]], t_srs_3763)
+	lims_3763 = (minimum(corners_3763[:, 1]), maximum(corners_3763[:, 1]),
+	             minimum(corners_3763[:, 2]), maximum(corners_3763[:, 2]))
+	I_native = wmsread(wms; layer=layer_n, region=lims_3763, size=(mMo*256, nMo*256))
+
+	# 2) Warp 3763 -> 3857 (Bing's intermediate mosaic) with the same `-te` / `-ts` as the Bing path.
+	xm, ym = geog2merc(lon_mm, lat_mm, pt_radius)
+	merc_srs::String = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=$pt_radius +b=$pt_radius +units=m +no_defs"
+	I_3857 = gdalwarp(I_native, ["-t_srs", merc_srs, "-r", "cubic", "-ts", string(nMo*256), string(mMo*256),
+	                              "-te", string(xm[1]), string(ym[1]), string(xm[2]), string(ym[2])])::GMTimage{UInt8, 3}
+
+	# 3) Final warp 3857 -> latlong (auto-size). Identical to Bing's last gdalwarp call.
+	I_geog = gdalwarp(I_3857, ["-t_srs", "+proj=latlong +datum=WGS84", "-r", "cubic"])::GMTimage{UInt8, 3}
+	return I_geog
+end
+
 # ---------------------------------------------------------------------------------------------------
 """
 	quadtree = XY2quadtree(x, y, zoom; quadkey::Matrix{Char}=['0' '1'; '2' '3']) -> String
@@ -587,7 +660,7 @@ Get information about a tile provider given its name and zoom level. The returne
 is only relevant for internal use and is an implementation detail not documented here.
 
 ### Arguments
-- `name`: Name of the tile provider. Currently available are "Bing" (the default), "Google", "OSM", "Esri", "Nimbo".
+- `name`: Name of the tile provider. Currently available are "Bing" (the default), "Google", "OSM", "Esri", "Nimbo", "ortosat23".
   Optionally, the `name` can be a tuple of two strings, where the first string is the provider name
   and the second string is the variant name (see the `variant` bellow).
 
@@ -605,6 +678,17 @@ is only relevant for internal use and is an implementation detail not documented
   - `Google`: variants => "Satellite", "Road", "Terrain", or "Hybrid".
   - `Esri`: variants => "World\\_Street\\_Map" (default), "Elevation/World\\_Hillshade", or "World\\_Imagery".
   - `Nimbo`: variants => "RGB" (default), "NIR", "NDVI", or "RADAR".
+  - `ortosat23` / `Ortosat23`: variants => "CorVerdadeira" (default) or "FalsaCor". Fetched via the DGT
+     WMS service (https://ortos.dgterritorio.gov.pt/wms/ortosat2023) — one WMS request covers the full
+     tile-matrix bbox. Always returned in WGS84 geographical coordinates. Spelling controls the pipeline:
+      * `"ortosat23"` (lowercase): single-warp approximate path. Output size slightly different from
+        same-zoom Bing/Google result (3763 vs 3857 aspect after warping to geog).
+      * `"Ortosat23"` (capital `O`): double-warp exact path. Two warps (3763 -> 3857 -> latlong)
+        replicate Bing's pipeline so output bbox/pixel count match Bing exactly.
+     For native EPSG:3763 metric coordinates or custom pixel resolution / `zoom` / `refine` factors,
+     call `wmsread` directly: first obtain the WMS handle with
+     `wms = wmsinfo("https://ortos.dgterritorio.gov.pt/wms/ortosat2023?service=wms&request=getcapabilities")`
+     and then `wmsread(wms; layer=..., region=..., pixelsize=...)`.
 """
 getprovider(name::Tuple{String,String}, zoom::Int; date::String="", key="") = getprovider(name[1], zoom; variant=name[2], date=date, key=key)
 function getprovider(name::StrSymb, zoom::Int; variant="", format="", ZYX::Bool=false, dir_code="", date::String="", key::String="")
@@ -649,11 +733,18 @@ function getprovider(name::StrSymb, zoom::Int; variant="", format="", ZYX::Bool=
 		variant = d * "_" * v * filesep * variant		# Need to carry also the dates (this will be used in cache name)
 		url = "https://prod-data.nimbo.earth/mapcache/tms/1.0.0/" * d * "_" * v * "@kermap/"
 		max_zoom = 13;	isZXY = true; code = "nimbo";	ext = "png";	sitekey = "?kermap_token=" * key
-	#elseif (_name == "ortosat")							# 
-		# variants: CorVerdadeira, falsacor
-		#layer = (variant == "falsacor") ? 1 : 2
-		#wms = wmsinfo("https://ortos.dgterritorio.gov.pt/wms/ortosat2023?service=wms&request=getcapabilities")
-		#url, isZXY, max_zoom, code = wms, true, 19, "ortosat"
+	elseif (_name == "ortosat23")						# Portuguese DGT WMS orthophotos 2023
+		# Tile XYZ is computed normally, but the actual image is fetched via WMS (wmsinfo/wmsread).
+		# variants: "CorVerdadeira" (default) or "FalsaCor".
+		# Provider name spelling controls the output pipeline:
+		#   "ortosat23" (lowercase) -> single-warp approximate path (output ~matches XYZ tile providers).
+		#   "Ortosat23" (capital O) -> double-warp exact path (output bbox/size exactly matches Bing).
+		(variant == "") && (variant = "CorVerdadeira")
+		vv = lowercase(variant)
+		(vv != "corverdadeira" && vv != "falsacor") &&
+			error("Ortosat23 only supports 'CorVerdadeira' (default) or 'FalsaCor' variants")
+		url = "https://ortos.dgterritorio.gov.pt/wms/ortosat2023?service=wms&request=getcapabilities"
+		max_zoom = 19;	isZXY = true;	code = (string(name)[1] == 'O') ? "Ortosat23" : "ortosat23"
 	elseif (_name == "mesh")							# For mesh we don't care the provider and max zoom is ilimitted
 		url, isZXY, max_zoom, code = string(name), true, 50, "unknown"
 	else