How to read a GeoTIFF after making a WMS request using Python

Question

I'm trying to read altitude GeoTIFFS from a WMS service. I know how to do this with BytesIO if the output format were a JPEG, but applying the same trick for rasterio does not seem to work. Does anyone have any suggestions?

url = 'http://geodata.nationaalgeoregister.nl/ahn3/wms?service=wms&version=1.3.0'
wms = WebMapService(url)

x1new= 51
x2new = 51.1
y1new = 5
y2new = 5.1

layer= 'ahn3_05m_dtm'



img = wms.getmap(layers = [layer], srs = 'EPSG:3857', bbox = [x1new,y1new,x2new,y2new] , size = (width,height), format= 'image/GeoTIFF')


r = rasterio.open(BytesIO(img.read()))
#this last step produces an error
r.read()

The last step hits me with an error

AttributeError: '_GeneratorContextManager' object has no attribute 'read'

score 2 · Answer 1 · answered Jul 17 '20 at 13:54

You probably already have found a solution, but for the record, this code would work. The problem was in the way tou opened your file.

url = 'http://geodata.nationaalgeoregister.nl/ahn3/wms?service=wms&version=1.3.0'
wms = WebMapService(url)

x1new= 100000
x2new = 100500
y1new = 450000
y2new = 450500

layer= 'ahn3_05m_dtm'


img = wms.getmap(layers = [layer], srs = 'EPSG:28992', bbox = [x1new,y1new,x2new,y2new] , size = (1920, 592), format= 'image/GeoTIFF')

with rasterio.open(BytesIO(img.read())) as r:
    thing = r.read()
    show(thing)

score 1 · Answer 2 · answered Jun 05 '20 at 09:03

You can use rasterio's MemoryFile class:

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

url = 'https://services.terrascope.be/wms/v2?'
wms = WebMapService(url)

x_min = 556945.9710290054
y_min = 6657998.9149440415
x_max = 575290.8578174476
y_max = 6663655.255037144

layer = 'CGS_S2_RADIOMETRY'

img = wms.getmap(
    layers = [layer],
    srs = 'EPSG:3857',
    bbox = (x_min, y_min, x_max, y_max),
    size = (1920, 592),
    format = 'image/png',
    time = '2020-06-01'
)

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

I didn't manage to reproduce your example though. Note that you seem to be using lat, long coordinates but querying the wms in EPSG:3857. Here is the output of the example above:

PS: for this kind of GIS-specific questions gis.stackexchange is probably more relevant.

How to read a GeoTIFF after making a WMS request using Python

2 Answers2