很久沒有寫博客了，確實(shí)有些懶。
今天抽空更新一次，下一次又不知道什么時間啦。。。

今天的主題是大家都已經(jīng)說了很多的白化，關(guān)于basemap的白化，最早氣象家園有個帖子講到，那后來用于cartopy的方法基本都是類似的邏輯和方法。
這幾天我也使用到了白化，但不同的是我在南海子圖中使用蘭伯特投影進(jìn)行白化，于是出現(xiàn)了非常奇怪的問題，下面我們慢慢講。

首先還是導(dǎo)入必要的庫

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文件也是從網(wǎng)上下載,現(xiàn)在已找不到出處,如有侵權(quán)敬請諒解.
同時也對原文件進(jìn)行了一定的修改,主要是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

設(shè)置下reload以便修改了函數(shù)隨時加載.

%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，發(fā)現(xiàn)在shaperead讀取時還存在問題，于是后面又重新轉(zhuǎn)換成字符串...
可以說這是一個不太友好的現(xiàn)象，盡管Python在推薦使用pathlib，但有些庫還沒有及時更新使用...

定義投影信息和測試數(shù)據(jù)

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)

定義函數(shù)將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

一切正常，圖片顯示的都是測試數(shù)據(jù)，那么下面開始白化大圖中的數(shù)據(jù)

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

正常,白化了大圖的數(shù)據(jù),下面白化南海子圖.

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

好了,我們發(fā)現(xiàn)了非常奇怪的現(xiàn)象,盡管白化了南海子圖中的數(shù)據(jù)，但是其他地區(qū)的數(shù)據(jù)也出現(xiàn)在了圖片中...

于是開始絞盡腦汁想解決辦法，最終想是否可以將白化的范圍限定在axis的extents內(nèi)呢?于是到matplotlib官網(wǎng)查看Path和Pathpatch的說明，終于找到了想要的結(jié)果，且看代碼.

#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
    '''
    當(dāng)region為None或不指定輸入時則為shp文件的全部屬性；
    可指定region參數(shù)為list或特殊的數(shù)值
    需要判斷所使用shp文件中的屬性(attr),然后修改shape_rec.record[*],此處的屬性與region對應(yīng)
    '''
    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匹配的唯一標(biāo)識符，record[]中必有一項(xiàng)是對應(yīng)的。
            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  重點(diǎn)看這里，增加了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

一切又正常了,多余不廢話了,有些代碼也是從網(wǎng)絡(luò)獲取,也看出一些更有的解決方法,但是懶得去改動了,直接使用把...

寫在最后，有個想法，可否針對cartpy的白化寫個開源庫，并把我們的地圖附上去，一方面推廣我們地圖邊界，另一方面也更加方便地使用白化。尋找志同道合的小伙伴一起抽時間做做把。