如何用圆形Patch裁剪瓦片影像？解决set_clip_path失效问题

问题描述 #

尝试用圆形Patch裁剪Cartopy加载的OSM卫星背景图，但调用osm_img.set_clip_path(patch)后无效果，原代码如下：

import io
from PIL import Image

import cartopy.crs as ccrs
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from urllib.request import urlopen, Request

def image_spoof(self, tile):
    '''this function reformats web requests from OSM for cartopy'''

    url = self._image_url(tile)                # get the url of the street map API
    req = Request(url)                         # start request
    req.add_header('User-agent','Anaconda 3')  # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read())            # get image
    fh.close()                                 # close url
    img = Image.open(im_data)                  # open image with PIL
    img = img.convert(self.desired_tile_form)  # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy

proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

# set extent
lon_min = -98.853627
lon_max = -98.752037
lat_min = 19.274685
lat_max = 19.376275

ax.set_extent((lon_min, lon_max, lat_min, lat_max), crs=proj)

cimgt.QuadtreeTiles.get_image = image_spoof # reformat web request for street map spoofing
osm_img = cimgt.QuadtreeTiles() # spoofed, downloaded street map

osm_img = ax.add_image(osm_img, 16) # add OSM with zoom specification

patch = patches.Circle((-98.802832, 19.32548), radius=0.03, transform=ax.transData)
ax.add_patch(patch)

# osm_img.set_clip_path(patch)  # 此代码无效

解决方案 #

原因分析 #

add_image返回的AxesImage对象对应分块加载的OSM瓦片，直接对其设置set_clip_path无法覆盖所有瓦片，因此裁剪效果不生效。

方法1：给坐标轴设置全局裁剪路径 #

直接将裁剪路径应用到整个坐标轴，所有绘图元素都会被圆形裁剪，适合只需要显示圆形区域内容的场景：

# ... 保留原代码中加载OSM和创建patch的部分 ...

patch = patches.Circle((-98.802832, 19.32548), radius=0.03, transform=ax.transData)
# 如果不需要显示圆形边框，设置透明填充和边框
patch.set_facecolor('none')
patch.set_edgecolor('none')
ax.add_patch(patch)

# 给坐标轴设置裁剪路径
ax.set_clip_path(patch)

plt.show()

方法2：手动合成裁剪后的OSM图像 #

如果需要仅裁剪OSM图层，保留坐标轴其他元素，可以先渲染OSM图像，再用PIL做圆形裁剪，最后重新绘制：

import io
from PIL import Image, ImageDraw

import cartopy.crs as ccrs
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from urllib.request import urlopen, Request

def image_spoof(self, tile):
    url = self._image_url(tile)
    req = Request(url)
    req.add_header('User-agent','Anaconda 3')
    fh = urlopen(req)
    im_data = io.BytesIO(fh.read())
    fh.close()
    img = Image.open(im_data)
    img = img.convert(self.desired_tile_form)
    return img, self.tileextent(tile), 'lower'

proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

lon_min = -98.853627
lon_max = -98.752037
lat_min = 19.274685
lat_max = 19.376275
ax.set_extent((lon_min, lon_max, lat_min, lat_max), crs=proj)

# 加载并绘制OSM图像
cimgt.QuadtreeTiles.get_image = image_spoof
osm_img = cimgt.QuadtreeTiles()
ax.add_image(osm_img, 16)

# 渲染当前画布到PIL图像
plt.savefig('temp_osm.png', bbox_inches='tight', pad_inches=0)
osm_pil = Image.open('temp_osm.png').convert('RGBA')

# 创建圆形蒙版
mask = Image.new('L', osm_pil.size, 0)
draw = ImageDraw.Draw(mask)
# 将经纬度坐标转换为图像像素坐标
center_x, center_y = ax.transData.transform((-98.802832, 19.32548))
# 计算半径对应的像素值（根据坐标轴范围和图像尺寸换算）
radius_pix = ax.transData.transform((-98.802832 + 0.03, 19.32548))[0] - center_x
draw.ellipse((center_x - radius_pix, center_y - radius_pix, 
              center_x + radius_pix, center_y + radius_pix), fill=255)

# 合成裁剪后的图像
cropped_osm = Image.composite(osm_pil, Image.new('RGBA', osm_pil.size, (255,255,255,0)), mask)

# 清空坐标轴，重新绘制裁剪后的图像
ax.clear()
ax.set_extent((lon_min, lon_max, lat_min, lat_max), crs=proj)
ax.imshow(cropped_osm, extent=(lon_min, lon_max, lat_min, lat_max), origin='upper', transform=proj)

# 可选：添加圆形边框
patch = patches.Circle((-98.802832, 19.32548), radius=0.03, transform=ax.transData, 
                       facecolor='none', edgecolor='red', linewidth=2)
ax.add_patch(patch)

plt.show()

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