[GIS] Reading a WMS/ raster map with rasterio

owspythonrasteriowms

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?

Best Answer

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 :)

enter image description here

Related Question