cartopy白化分享_1

坎坷

本人也是从论坛中和网上扒的python白化教程，应该用的都是一个方法，这里不再赘述。

示例脚本：
# -*- encoding: utf-8 -*-
import maskout_new
import numpy as np
import lambert_ticks
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER

# 站点数据
stid_lon = [89.2, 82.95, 75.98, 79.93, 85.55, 88.17, 93.52, 88.08, 83, 84.66, 81.33, 87.62]
stid_lat = [42.93, 41.72, 39.47, 37.13, 38.15, 39.03, 42.82, 47.73, 46.73, 44.43, 43.95, 43.78]
stid_data = np.random.rand(len(stid_lon)) * 500

# 站点数据格点化
nlon = np.linspace(70, 100, 90)
nlat = np.linspace(30, 55, 75)
nlon, nlat = np.meshgrid(nlon, nlat)
func = Rbf(stid_lon, stid_lat, stid_data, function='linear')
stid_pre = func(nlon, nlat)

# 画图
proj = ccrs.LambertConformal(central_latitude=90, central_longitude=84.6)
fig, ax = plt.subplots(figsize=(15, 15), subplot_kw=dict(projection=proj))  # 建立页面

# 白化
shpfile = "/home/frg/Python3_Data/Data/shp_file/China_shapefiles/Province.shp"  # 读取目标shp
area_str = ["新疆维吾尔自治区"]  # 目标区域
ax.set_extent([74, 95.2, 33.8, 49], ccrs.PlateCarree())  # 设置经纬度范围
con = ax.contourf(nlon, nlat, stid_pre, transform=ccrs.PlateCarree())  # 等值线图
con_mask = maskout_new.shp2clip(con, ax, shpfile=shpfile, area_str=area_str, proj=proj)
# 站点标记
ax.scatter(stid_lon, stid_lat, c='w', s=25, marker='o', transform=ccrs.PlateCarree())

# 坐标与经纬网格（兰伯特投影）
fig.canvas.draw()
xticks = list(range(-180, 180, 4))
yticks = list(range(-90, 90, 3))
# ax.gridlines(xlocs=xticks, ylocs=yticks, linewidth=1.2, linestyle='--')  # 经纬线网格
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_ticks.lambert_xticks(ax, xticks)
lambert_ticks.lambert_yticks(ax, yticks)

# 出图
# plt.show()
plt.savefig('/home/frg/Python3_Data/Figure/map_1.png')

这个示例脚本也属于缝合怪--是各大网站论坛中的帖子扒下来的。但其中引用的maskout_new为笔者自己改编，其中添加了一些实用功能：
maskout_new.py内容如下：

# -*- encoding: utf-8 -*-
import shapefile
import cartopy.crs as ccrs
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from cartopy.io.shapereader import Reader


def find_region(shp_path, area_str):
    file = shapefile.Reader(shp_path)  # 读取
    opt = 0
    for k in range(len(file.records()[0])):
        i = 0
        for shape_rec in file.shapeRecords():
            if isinstance(shape_rec.record[k], bytes):
                if str(shape_rec.record[k].decode("gbk")) == area_str:
                    return k, i, shape_rec.record[k]
            i = i + 1
    if opt == 0:
        print("无此区域")
        return 0, 0


def shp2clip(originfig, ax, shpfile, area_str, proj=None):
    region = []
    for i in range(len(area_str)):
        region_k, region_i, region_str = find_region(shpfile, area_str)
        ax.add_geometries(list(Reader(shpfile).geometries())[region_i:region_i + 1], ccrs.PlateCarree(),
                          facecolor='none', edgecolor='w', linewidth=0.5)
        region.append(region_str)
    sf = shapefile.Reader(shpfile)
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        if shape_rec.record[int(region_k)] in 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, 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 - 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



maskout_new可实现以下功能：
1.实现汉字输入查找区域。

1. area_str = ["新疆维吾尔自治区"]  # 目标区域

复制代码

你也可以输入

1. area_str = ["新疆维吾尔自治区","内蒙古自治区"]  # 目标区域

复制代码

即可简便的白化想要区域
其中原理就是通过“gbk”转码，将shp文件中关键字转换为汉字再匹配相应区域

1. if str(shape_rec.record[k].decode("gbk")) == area_str:

复制代码

2.尽量实现读取所有来源的shp文件。由于网上的shp文件繁多，每个shp文件属性不一样，现在看来通过匹配汉字是一个通用方法。笔者已经测试通过micaps目录下的shp文件（如Province.shp）、国家基础地理数据（如bou2_4p.shp）等都可以顺利读取，如省级shp文件可以通过area_str输入如“四川省”等白化所需区域。更加细致的shp文件如市级、县级，也可以通过area_str输入某某市、某某县白化所需区域，但要注意是，如果白化区域过多，脚本运行时间可能较长。

3.自动匹配region_k，这个参数可以理解为查找的区域的关键字在shp文件属性中的第几列（大致可以这样理解）。在之前我扒的帖子中这个参数是定死的，如3或4又或7。但问题是读其他与作者不相同的shp文件时，这个参数肯定是不同的，但小白又不清楚问题在哪儿。所以笔者在匹配区域时自动返回region_k，并调用shp2clip函数（非原创）对图像进行裁剪。

最后。。。
使用说明：
更改以下参数，即可快捷的白化所需区域

1. shpfile = "/home/frg/Python3_Data/Data/shp_file/China_shapefiles/Province.shp"  # 读取目标shp
   
2. area_str = ["新疆维吾尔自治区"]  # 目标区域

复制代码

zzw9701

为什么我在运行的时候提示：

  File "D:\chengxu\Python\EXP\read_wrf.py", line 56, in <module>
    maskout_new.shp2clip(rh,ax1,provinces_shp,["广东省","广西省","福建省"],proj=proj)

  File "D:\chengxu\Python\EXP\maskout_new.py", line 26, in shp2clip
    region_k, region_i, region_str = find_region(shpfile, area_str)

ValueError: not enough values to unpack (expected 3, got 2)
就是在您的这行代码中要输出三个值，可是只输入了2个值
region_k, region_i, region_str = find_region(shpfile, area_str)

zhongyon

Asofhknok 发表于 2020-8-18 22:57
好人一生平安，问一个其他的，请问楼主知不知道怎么根据shp判断某格点在不在范围内呢？

import shapefile
import shapely.geometry as geometry

lon = 104
lat = 35

shp_file = shapefile.Reader('shp/aa.shp', encoding='gbk')
shp_file = shp_file.shapeRecords()[0].shape

if geometry.Point(lon, lat).within(geometry.shape(shp_file)):
    print('点在shp里面')

Asofhknok

好人一生平安，问一个其他的，请问楼主知不知道怎么根据shp判断某格点在不在范围内呢？

紫依

谢谢分享，正在学习

向南V5

zhongyon 发表于 2020-8-19 08:54
import shapefile
import shapely.geometry as geometry

这个用起来感觉效率有点低，大佬有没有其他的

华仔_beyond

感谢分享！
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xba in position 0: invalid start byte
大佬，请问一下这个是怎么回事呢？

zhongyon

华仔_beyond 发表于 2020-8-19 13:36
感谢分享！
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xba in position 0: invalid start by ...

你shp文件是gbk编码的吧

华仔_beyond

zhongyon 发表于 2020-8-19 16:27
你shp文件是gbk编码的吧

具体我也不清楚，我用的是'bou2_4p.shp'文件，用arcgis10.2打开是属性表中文，没有乱码。。我该怎么操作

Asofhknok

zhongyon 发表于 2020-8-19 08:54
import shapefile
import shapely.geometry as geometry

谢谢大佬！

guozhaodong221

请问楼主lambert_ticks这个库怎么安装，我conda和pip试过都不行