basemap 画图按照 shape 文件裁切图片，经纬度转为web墨卡托

相关文章推荐

腼腆的土豆 · WRF学习——使用CMIP6数据驱动WRF/ ...· 1 周前 ·

谈吐大方的花卷 · 全球法规网－汤加· 2 周前 ·

爱听歌的铁板烧 · 过程装备与控制工程_百度百科· 3 月前 ·

谦虚好学的太阳 · VG官宣打野离队，曾2次单杀王思聪的男人，如 ...· 10 月前 ·

豪情万千的台灯 · 产品经理为什么一定要学paddlepaddl ...· 1 年前 ·

豁达的圣诞树 · 当宠物医生需要考什么证吗？· 1 年前 ·

   vis_clevs = [0, 0.2, 0.5, 1, 2, 3, 5, 10, 20, 30, 9999],
   vis_colors = ['#642B00', '#9A00FF', '#FE0300', '#F95C00', '#FDB157', '#FFFF00', '#76FD2C', '#94E0EC', '#C7ECFF', '#FFFFFF']


    fig = plt.figure(figsize=(16, 8))
    ax = fig.add_axes([-0.00001, 0, 1.00001, 1])
    sf = shapefile.Reader(shppath + '.shp') # encoding='utf-8'  默认 utf-8  也可以设置别的
    m = Basemap(llcrnrlon=left_lon, urcrnrlon=right_lon,
                llcrnrlat=left_lat, urcrnrlat=right_lat, resolution='l', ax=ax)
　　# 添加自定义 色卡 和值
    cmap_mosaic = mpl.colors.ListedColormap(colors)
    norm = mpl.colors.BoundaryNorm(clevss, ncolors=cmap_mosaic.N, clip=True)
    im1 = m.contourf(lons, lats, data, levels=clevss, norm=norm, cmap=cmap_mosaic)
    clip, vertices = clip_region(sf, ax)
    for cour in im1.collections:
        cour.set_clip_path(clip)

def clip_region(sf, ax):
    把地图按照shape文件裁切
    :param sf:
    :param ax:
    :return:
    for shape_rec in sf.shapeRecords():
        vertices = []
        codes = []
        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]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [matplotlib.path.Path.MOVETO]
            codes += [matplotlib.path.Path.LINETO] * (prt[i + 1] - prt[i] - 2)
            codes += [matplotlib.path.Path.CLOSEPOLY]
        clip = matplotlib.path.Path(vertices, codes)
        clip = matplotlib.patches.PathPatch(clip, transform=ax.transData)
    return clip, vertices

def wgs84toWebMercator(lon, lat):
    将经纬度转为web墨卡托下的坐标,转换公式如下
    :param lon: 经度,需传入numpy数组,范围为[-180,180]
    :param lat: 纬度,需传入numpy数组,范围为[-85.05，85.05]
    :return:x,横坐标,对应于经度;y,纵坐标,对应于纬度
    x = lon * 20037508.342789 / 180
    y = np.log(np.tan((90 + lat) * np.pi / 360)) / (np.pi / 180)
    y = y * 20037508.34789 / 180
    return x, y
def EqlLatLonToWebMerc(sPicName, nIsLevel=0):
        将等经纬度图转为web墨卡托投影图
        :param sPicName: 传入等经纬度图像全路径名称
        :param nIsLevel: 是否为EC分层绘图
    image = Image.open(sPicName)
    data = np.array(Image.open(sPicName))
    lon, lat = np.meshgrid(np.linspace(left_lon, right_lon, image.size[0]),
                           np.linspace(right_lat, left_lat, image.size[1]), )
    color_r = data[:, :, 0]
    color_g = data[:, :, 1]
    color_b = data[:, :, 2]
    color_a = data[:, :, 3]
    lon_0, lon_1, lat_0, lat_1 = 60, 150, 0, 60  # 经纬度范围
    if nIsLevel == 1:
        resolution = 8000  # 设置分辨率
    else:
        resolution = 4000
    all_x, all_y = wgs84toWebMercator(
        np.array([lon_0, lon_1]), np.array([lat_0, lat_1]))
    all_col = (all_x[1] - all_x[0]) / resolution
    all_row = (all_y[1] - all_y[0]) / resolution
    x, y = wgs84toWebMercator(lon, lat)
    target_col = ((x - all_x[0]) / resolution)[0]
    target_row = ((all_y[1] - y) / resolution)[:, 0]
    f = insert_data2d(target_col, target_row, color_r, )
    res_r = f(np.arange(int(all_col)), np.arange(int(all_row)))
    f = insert_data2d(target_col, target_row, color_g, )
    res_g = f(np.arange(int(all_col)), np.arange(int(all_row)))
    f = insert_data2d(target_col, target_row, color_b, )
    res_b = f(np.arange(int(all_col)), np.arange(int(all_row)))
    f = insert_data2d(target_col, target_row, color_a, )
    res_a = f(np.arange(int(all_col)), np.arange(int(all_row)))
    Image.fromarray(np.dstack((res_r, res_g, res_b, res_a)
                              ).astype('uint8')).save(sPicName)