Reading a WMS/ raster map with rasterio

Question

Can I read a raster map which has been constructed in mapserver as WMS with rasterio?

I tried :

from rasterio.io import MemoryFile

def print_metadata(dataset):
     print(dataset.profile)


url ='http://localhost:8080/?map=/maps/ivm3.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=tr_ivm_gray&STYLES=&SRS=epsg:3857&BBOX=2857613.741389,4274927.875099,4989229.633477,5176940.449967&WIDTH=1160&HEIGHT=540&FORMAT=image/tiff'

tif_bytes = open(url,'rb').read()

with MemoryFile(tif_bytes) as memfile:
     with memfile.open() as dataset:
         print_metadata(dataset)

and it gives :

FileNotFoundError: [Errno 2] No such file or directory: 'http://localhost:8080/?map=/maps/ivm3.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=tr_ivm_gray&STYLES=&SRS=epsg:3857&BBOX=2857613.741389,4274927.875099,4989229.633477,5176940.449967&WIDTH=1160&HEIGHT=540&FORMAT=image/tiff'

another code :

raster=rasterio.open(url)
print (type(raster))

print ( raster.meta )
print ( raster.count )

array=raster.read()

this gives:

<class 'rasterio.io.DatasetReader'>
{'driver': 'WMS', 'dtype': 'uint8', 'nodata': None, 'width': 1073741824, 'height': 454363579, 'count': 3, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.001985222000710666, 0.0, 2857613.741389,
       0.0, -0.0019852220040462356, 5176940.449967)}
3
Traceback (most recent call last):
  File "rio1.py", line 20, in <module>
    array=raster.read()
  File "rasterio/_io.pyx", line 322, in rasterio._io.DatasetReaderBase.read
MemoryError

raster.read() tries to read too much data. What can I do here?

Answer 1

In case someone else runs into the same difficulty.

The MemoryFile class is indeed the way to go with rasterio.

In the first example you use MemoryFile correctly, but you're trying to request data from an url using python's open() built-in function, which is not possible as far as I know. The second example seems correct, but judging from the memory error you get and as indicated by @user2856, you're probably just trying to read too much data.

Here are a couple of working examples.

1) Reproducing the first example in the question using urllib and rasterio.MemoryFile:

from rasterio import MemoryFile
from rasterio.plot import show
from urllib.request import urlopen

url ='https://services.terrascope.be/wms/v2?service=WMS&version=1.3.0&request=GetMap&layers=CGS_S2_RADIOMETRY&format=image/png&time=2020-06-01&width=1920&height=592&bbox=556945.9710290054,6657998.9149440415,575290.8578174476,6663655.255037144&styles=&srs=EPSG:3857'

tif_bytes = urlopen(url).read()

with MemoryFile(tif_bytes) as memfile:
     with memfile.open() as dataset:
            print(dataset.profile)
            show(dataset)

2) Alternatively you can use owslib, which comes with a bunch of functionalities to query OGC web services:

from owslib.wms import WebMapService
from rasterio import MemoryFile
from rasterio.plot import show

wms_url = 'https://services.terrascope.be/wms/v2?'

wms = WebMapService(wms_url)

request = wms.getmap(
    layers=['CGS_S2_RADIOMETRY'],
    srs='EPSG:3857',
    format='image/png',
    bbox=(556945.9710290054,6657998.9149440415,575290.8578174476,6663655.255037144),
    size=(1920,592),
    time='2020-06-01'
)

with MemoryFile(request) as memfile:
     with memfile.open() as dataset:
            show(dataset)

3) Finally you can also directly read the data with rasterio.open(), as in the second example in the question:

import rasterio
from rasterio.plot import show

url ='https://services.terrascope.be/wms/v2?service=WMS&version=1.3.0&request=GetMap&layers=CGS_S2_RADIOMETRY&format=image/png&time=2020-06-01&width=1920&height=592&bbox=556945.9710290054,6657998.9149440415,575290.8578174476,6663655.255037144&styles=&srs=EPSG:3857'

raster = rasterio.open(url)

print(raster.meta)
show(raster)

And here is the result :)

Answer 2

I had the same problem, and I solved it mostly thanks to Read from TIF overviews using rasterio.

My XML file looks like this:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://a.tiles.wmflabs.org/bw-mapnik/${z}/${x}/${y}.png</ServerUrl>
        <SRS>EPSG:3857</SRS>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache>
        <Path>../.gdalwmscache</Path>
        <Extension>.png</Extension>
    </Cache>
</GDAL_WMS>

But any GDAL_WMS file should do. Note the TileLevel parameter, this sets the highest zoom level you want to go to. So even if this tileset supported zoom level 20, GDAL/rasterio will only see as far as 18.

The basic Python code I have is:

TMS_dataset = rasterio.open('templates/gdal_tms.xml')
# Calculate the pixel positions of the desired bounding box at the highest zoom level
# as specified in the XML file.
bl = TMS_dataset.index(west, south, op=math.floor)
tr = TMS_dataset.index(east, north, op=math.ceil)
# image_size is a tuple (h, w, num_bands)
output_dataset = np.empty(shape=image_size, dtype=TMS_dataset.profile['dtype'])

# Read each band
TMS_dataset.read(1, out=output_dataset[:, :, 0], window=((tr[0], bl[0]), (bl[1], tr[1])))
TMS_dataset.read(2, out=output_dataset[:, :, 1], window=((tr[0], bl[0]), (bl[1], tr[1])))
TMS_dataset.read(3, out=output_dataset[:, :, 2], window=((tr[0], bl[0]), (bl[1], tr[1])))

# Create an output image dataset
output_image = rasterio.open('flange3.png', 'w', driver='png', width=image_size[1], height=image_size[0], count=3,
                             dtype=output_dataset.dtype)
# Write each band
output_image.write(output_dataset[:, :, 0], 1)
output_image.write(output_dataset[:, :, 1], 2)
output_image.write(output_dataset[:, :, 2], 3)

output_image.close()

This particular code doesn't bother with georeferencing the output because I didn't need it for this application, but the rasterio docs tell you how to achieve that if you need it. It will emit a NotGeoReferencedwarning and will complain that the output file doesn't exist, but they are harmless as far as I can tell.

Note that it doesn't read in every tile at every level, GDAL/rasterio is clever enough to read just the tiles that meet the desired resolution. In addition, adding the Cache option to the XML caches tiles locally. This is only valid for TMS-like datasets, but is absolutely perfect for my use-case.

In addition, despite the docs saying that you can't resample with "decimated reads", you can add resampling=Resampling.cubic to the dataset read() function and it does the right thing. Although this might only apply to this use-case.