Search code examples
c++webgdalimaging

GDAL read several pictures from wmts server using same opened connection

I use C++ code to read pictures from WMTS server using DGAL.

First I initialize GDAL once:

...
OGRRegisterAll();
etc.

But new connection is opened every time I want to read new image (different urls):

gdalDataset = GDALOpen(my_url, GA_ReadOnly);

URL example: https://sampleserver6.arcgisonline.com/arcgis/rest/services/Toronto/ImageServer/tile/12/1495/1145

Unfortunately I didn't find a way to read multiply images by same connection.

Is there such option in GDAL or in WMTS? Are there other ways to improve timing (I read thousands of images)?

Solution

  • While GDAL can read PNG files, it doesn't add much since those lack any geographical metadata.

    You probably want to interact with the WMS server instead, not the images directly. You can for example run gdalinfo on the main url to see the subdatasets:
    gdalinfo https://sampleserver6.arcgisonline.com/arcgis/services/Toronto/ImageServer/WMSServer?request=GetCapabilities&service=WMS

    The first layer seems to have an issue, I'm not sure, but the other ones seem to behave fine.

    I hope you don't mind me using some Python code, but the c++ api should be similar. Or you could try using the command-line utilities first (gdal_translate), to get familiar with the service.

    See the WMS driver for more information and examples:
    https://gdal.org/drivers/raster/wms.html

    You can for example retrieve a subset and store it with:

    from osgeo import gdal

url = r"WMS:https://sampleserver6.arcgisonline.com:443/arcgis/services/Toronto/ImageServer/WMSServer?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=Toronto%3ANone&SRS=EPSG:4326&BBOX=-79.454856,43.582524,-79.312167,43.711781"

bbox = [-79.35, 43.64, -79.32, 43.61]
filename = r'D:\Temp\toronto_subset.tif'

ds = gdal.Translate(filename, url, xRes=0.0001, yRes=0.0001, projWin=bbox)
ds = None

    Which looks like:

    import numpy as np
import matplotlib.pyplot as plt

ds = gdal.OpenEx(filename)
img = ds.ReadAsArray()
ds = None

mpl_extent = [bbox[i] for i in [0,2,3,1]]

fig, ax = plt.subplots(figsize=(5,5), facecolor="w")
ax.imshow(np.moveaxis(img, 0, -1), extent=mpl_extent)

    enter image description here

    Note that the data in native resolution for these type of services is often ridiculously large, so usually you want to specify a subset and/or limited resolution as the output.