Python从WMS服务器下载的TIFF文件数值与原文件不符问题排查

使用Python通过WMS服务下载GeoTIFF文件，与网站直接下载的原文件对比，二者尺寸、目标坐标的行列号完全一致，但对应位置的数值却不匹配：

下载代码 #

import requests
import numpy as np
import rasterio
from rasterio.transform import from_origin
from rasterio.transform import rowcol
from rasterio.merge import merge
from rasterio.plot import show
from pyproj import Transformer
import certifi
import ssl
import matplotlib.pyplot as plt

Layer = "warme_nachten_huidig"
minx=13588.806396509579
miny=306893.81943950325
maxx=277788.8063965096
maxy=612693.8194395032
map_width = "2642"
map_height = "3058"
bbox = (minx, miny, maxx, maxy)

def map_downloaden(bbox, width, height):
    ssl.SSLContext.verify_mode = ssl.VerifyMode.CERT_OPTIONAL

    wms_url = "https://apps.geodan.nl/public/data/org/gws/YWFMLMWERURF/kea_public/ows?"
    
    params = {
        "service": "WMS",
        "version": "1.3.0",
        "request": "GetMap",
        "layers": f"{Layer}",
        "styles": "",
        "crs": "EPSG:28992",
        "bbox": ','.join(map(str, bbox)),
        "width": width,
        "height": height,
        "format": "image/geotiff"
    }

    response = requests.get(wms_url, params=params, stream=True, verify=certifi.where())

    if response.status_code == 200:
        with open(file_path_python, 'wb') as f:
            f.write(response.content)
        print(f"Map successfully downloaded to {file_path_python}")
    else:
        print(f"Failed to download map. Status code: {response.status_code}")

map_downloaden(bbox, map_width, map_height)

数值对比代码 #

def get_coordinate_value(coord, file_path, versie):
    with rasterio.open(file_path) as dataset:
        transformer = Transformer.from_crs("EPSG:4326", dataset.crs, always_xy=True)
        projected_coord = transformer.transform(coord[0], coord[1])
        row, col = rowcol(dataset.transform, projected_coord[0], projected_coord[1])

        if 0 <= row < dataset.height and 0 <= col < dataset.width:
            data_value = dataset.read(1)[row, col]
            print(f"---- {versie}-------- Height: {dataset.height} Width: {dataset.width} ---- row: {row} col: {col} ---- COORDINATES: {coord} ---- VALUE: {data_value} ")
        else:
            print("Coordinates are out of the raster's bounds.")

coord = ( 5.975079577360298, 51.6954090242412)
get_coordinate_value(coord, file_path_python, "Python---")
get_coordinate_value(coord, file_path_original, "Original")

对比结果 #

---- Python----------- 
Height: 3058 Width: 2642 ---- row: 2006 col: 1820 ---- COORDINATES: (5.975079577360298, 51.6954090242412) ---- VALUE: 3 

---- Original--------- 
Height: 3058 Width: 2642 ---- row: 2006 col: 1820 ---- COORDINATES: (5.975079577360298, 51.6954090242412) ---- VALUE: 0.6150805354118347 

原因分析 #

• WMS默认返回渲染后的数据：WMS的GetMap请求默认会应用图层预设的渲染样式，把原始浮点型的连续数值转换为离散的分级整数（比如示例中的3是分级后的类别编号），而网站直接下载的是原始浮点数据。
• 数据类型量化：WMS返回的GeoTIFF可能被转换为8位/16位整数格式以适配图像渲染，而原文件保留了原始的浮点型数据，导致数值映射偏差。
• 请求参数未指定原始数据格式：当前请求仅指定format=image/geotiff，未明确要求获取原始未渲染的栅格数据，服务默认返回经过处理的图像数据。

解决办法 #

1. 修改WMS请求参数，获取原始数据
   调整params中的参数，明确请求原始浮点数据：

   params = {
       # 原有参数不变
       "format": "image/geotiff;subtype=GeoTIFF",
       "format_options": "FLOAT32=TRUE",
       "styles": "raster"  # 部分服务需要指定无渲染的样式，或留空
   }
   

   注：具体参数需参考该WMS服务的GetCapabilities文档（访问wms_url后添加request=GetCapabilities查看），确认支持的原始数据格式和参数。

2. 检查并确认数据类型
   先查看两个文件的数据类型差异，验证是否存在量化问题：

   with rasterio.open(file_path_python) as ds:
       print("Python下载文件数据类型:", ds.dtypes[0])
       print("Python下载文件数值范围:", ds.read(1).min(), ds.read(1).max())
   with rasterio.open(file_path_original) as ds:
       print("原文件数据类型:", ds.dtypes[0])
       print("原文件数值范围:", ds.read(1).min(), ds.read(1).max())
   

3. 改用WCS服务获取原始栅格
   如果该服务提供WCS（Web Coverage Service）接口，建议直接使用WCS的GetCoverage请求获取原始栅格数据，WCS专门用于获取未经过渲染的原始栅格，避免WMS的渲染处理。

内容的提问来源于stack exchange，提问作者user7342931