很久没有写博客了，确实有些懒。
今天抽空更新一次，下一次又不知道什么时间啦。。。
今天的主题是大家都已经说了很多的白化，关于basemap的白化，最早气象家园有个帖子讲到，那后来用于cartopy的方法基本都是类似的逻辑和方法。
这几天我也使用到了白化，但不同的是我在南海子图中使用兰伯特投影进行白化，于是出现了非常奇怪的问题，下面我们慢慢讲。
首先还是导入必要的库
import os
import pathlib
import datetime
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
import cartopy.io.shapereader as shpreader
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import maskout
%matplotlib inline
maskout文件也是从网上下载,现在已找不到出处,如有侵权敬请谅解.
同时也对原文件进行了一定的修改,主要是region可以为None,下面贴出代码。
#coding=utf-8
import shapefile
import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def shp2clip(originfig, ax, shpfile, region, proj=None):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if region is None:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i + 1]):
if proj:
vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.Geodetic()))
else:
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
codes += [Path.CLOSEPOLY]
continue
if shape_rec.record[3] == region:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i + 1]):
if proj:
vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.Geodetic()))
else:
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)
for contour in originfig.collections:
contour.set_clip_path(clip)
return clip
设置下reload以便修改了函数随时加载.
%load_ext autoreload
%autoreload 2
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
定义shapefile的路径，这些shp文件都是从meteoinfo的map下拷贝的。
homedir = os.getcwd()
tophome = pathlib.Path(homedir).parent
shpdir = tophome / 'china_shp'
cnpdir = shpdir / 'cn_province.shp'
ninedir = shpdir / 'cn_border.shp'
(cnpdir.exists(), ninedir.exists())
(True, True)
使用pathlib的Path，发现在shaperead读取时还存在问题，于是后面又重新转换成字符串...
可以说这是一个不太友好的现象，尽管Python在推荐使用pathlib，但有些库还没有及时更新使用...
定义投影信息和测试数据
lcc_proj = ccrs.LambertConformal(central_longitude=102, central_latitude=35)
xlon = np.arange(70,141,5)
ylat = np.arange(0,61,5)
glon, glat = np.meshgrid(xlon, ylat)
data = np.random.randint(0,10,glon.shape)
定义函数将shapefile加载到axis中
def add_shp(ax, shpdir, edgecolor):
reader = shpreader.Reader(shpdir)
shpdata = cfeature.ShapelyFeature(reader.geometries(), crs=ccrs.PlateCarree(), edgecolor=edgecolor, facecolor='none')
ax.add_feature(shpdata)
第一个测试图
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection=lcc_proj, position=[0.04,0.04,0.95,0.95])
ax.set_extent([70,140,10,55])
ax2 = fig.add_axes([0.1,0.2,0.3,0.3],projection=lcc_proj)
ax2.set_extent([106, 121, 5, 23], crs=ccrs.PlateCarree())
add_shp(ax, str(cnpdir), 'k')
add_shp(ax, str(ninedir), 'k')
add_shp(ax2, str(cnpdir), 'k')
add_shp(ax2, str(ninedir), 'k')
cs = ax.contourf(glon, glat, data, transform=ccrs.PlateCarree())
cs2 = ax2.contourf(glon, glat, data, transform=ccrs.PlateCarree())
# clip = maskout.shp2clip(cs, ax, str(cnpdir), None, lcc_proj)
# clip2 = maskout.shp2clip(cs2, ax2, str(cnpdir), None, lcc_proj)
output_12_0.png
一切正常，图片显示的都是测试数据，那么下面开始白化大图中的数据
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection=lcc_proj, position=[0.04,0.04,0.95,0.95])
ax.set_extent([70,140,10,55])
ax2 = fig.add_axes([0.1,0.2,0.3,0.3],projection=lcc_proj)
ax2.set_extent([106, 121, 5, 23], crs=ccrs.PlateCarree())
add_shp(ax, str(cnpdir), 'k')
add_shp(ax, str(ninedir), 'k')
add_shp(ax2, str(cnpdir), 'k')
add_shp(ax2, str(ninedir), 'k')
cs = ax.contourf(glon, glat, data, transform=ccrs.PlateCarree())
cs2 = ax2.contourf(glon, glat, data, transform=ccrs.PlateCarree())
clip = maskout.shp2clip(cs, ax, str(cnpdir), None, lcc_proj)
# clip2 = maskout.shp2clip(cs2, ax2, str(cnpdir), None, lcc_proj)
output_14_0.png
正常,白化了大图的数据,下面白化南海子图.
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection=lcc_proj, position=[0.04,0.04,0.95,0.95])
ax.set_extent([70,140,10,55])
ax2 = fig.add_axes([0.1,0.2,0.3,0.3],projection=lcc_proj)
ax2.set_extent([106, 121, 5, 23], crs=ccrs.PlateCarree())
add_shp(ax, str(cnpdir), 'k')
add_shp(ax, str(ninedir), 'k')
add_shp(ax2, str(cnpdir), 'k')
add_shp(ax2, str(ninedir), 'k')
cs = ax.contourf(glon, glat, data, transform=ccrs.PlateCarree())
cs2 = ax2.contourf(glon, glat, data, transform=ccrs.PlateCarree())
clip = maskout.shp2clip(cs, ax, str(cnpdir), None, lcc_proj)
clip2 = maskout.shp2clip(cs2, ax2, str(cnpdir), None, lcc_proj)
output_16_0.png
好了,我们发现了非常奇怪的现象,尽管白化了南海子图中的数据，但是其他地区的数据也出现在了图片中...
于是开始绞尽脑汁想解决办法，最终想是否可以将白化的范围限定在axis的extents内呢?于是到matplotlib官网查看Path和Pathpatch的说明，终于找到了想要的结果，且看代码.
#coding=utf-8
###################################################################################################################################
#####This module enables you to maskout the unneccessary data outside the interest region on a matplotlib-plotted output instance
####################in an effecient way,You can use this script for free ########################################################
#####################################################################################################################################
#####USAGE: INPUT include 'originfig':the matplotlib instance##
# 'ax': the Axes instance
# 'shapefile': the shape file used for generating a basemap A
# 'region':the name of a region of on the basemap A,outside the region the data is to be maskout
# OUTPUT is 'clip' :the the masked-out or clipped matplotlib instance.
from typing import List
import shapefile
import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def shp2clip(originfig,ax,shpfile,region=None, proj=None):#,region
'''
当region为None或不指定输入时则为shp文件的全部属性；
可指定region参数为list或特殊的数值
需要判断所使用shp文件中的属性(attr),然后修改shape_rec.record[*],此处的属性与region对应
'''
if isinstance(region, List):
region_type = 'in'
else:
region_type = 'target'
if region is None:
region_type = None
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if region_type is None:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
if proj:
vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.Geodetic()))
else:
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
continue
if region_type == 'in' and shape_rec.record[6] in region:
print(f'multi polygons : {region}')
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
if proj:
vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.PlateCarree()))
else:
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
elif region_type == 'target' and shape_rec.record[6] == region: ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
if proj:
vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.Geodetic()))
else:
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
# clip to ax.get_extent bbox 重点看这里，增加了Path首先clip到地图边界
bbox = ax.get_extent()
clip = clip.clip_to_bbox([bbox[0], bbox[2], bbox[1], bbox[3]], inside=True)
clip = PathPatch(clip, transform=ax.transData)
for contour in originfig.collections:
contour.set_clip_path(clip)
return clip
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection=lcc_proj, position=[0.04,0.04,0.95,0.95])
ax.set_extent([70,140,10,55])
ax2 = fig.add_axes([0.1,0.2,0.3,0.3],projection=lcc_proj)
ax2.set_extent([106, 121, 5, 23], crs=ccrs.PlateCarree())
add_shp(ax, str(cnpdir), 'k')
add_shp(ax, str(ninedir), 'k')
add_shp(ax2, str(cnpdir), 'k')
add_shp(ax2, str(ninedir), 'k')
cs = ax.contourf(glon, glat, data, transform=ccrs.PlateCarree())
cs2 = ax2.contourf(glon, glat, data, transform=ccrs.PlateCarree())
clip = shp2clip(cs, ax, str(cnpdir), None, lcc_proj)
clip2 = shp2clip(cs2, ax2, str(cnpdir), None, lcc_proj)
output_19_0.png
一切又正常了,多余不废话了,有些代码也是从网络获取,也看出一些更有的解决方法,但是懒得去改动了,直接使用把...
写在最后，有个想法，可否针对cartpy的白化写个开源库，并把我们的地图附上去，一方面推广我们地图边界，另一方面也更加方便地使用白化。寻找志同道合的小伙伴一起抽时间做做把。
