basemap 画图按照 shape 文件裁切图片, 经纬度转为web墨卡托
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)
参考网址:(比上面的还好)
https://blog.csdn.net/weixin_46604505/article/details/118770298