Python

2023-11-30 01:44| 来源: 网络整理| 查看: 265

在学习画图的过程中，看了许多大佬的绘图代码收益匪浅。在巨人的肩膀上继续前进，分享这一次的画图。多数没有注释，原理可能需要额外找别的帖子进行查阅。

再次之前，anaconda安装cartopy包也遇到了不少困难，我的解决方案是：装好对应Python版本的四个包：pyshp, Pillow, pyproj, Shapely。另外安装xarray的时候记得安装netcdf4。对应的安装网站在评论区进行分享。

本次画图笔记分成三个部分：

1.文件读取（采用的是WRF输出的wrfout文件）。

2.利用cartopy-读取shp文件，在前面几篇大佬的基础上动用掩模，以及投影兰伯特投影。这里用的是封装函数，我直接写在总代码里。

3.绘制场图，在colorbar上分成几种不同的形式。

读取包

import cartopy.crs as ccrs from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter from cartopy.io.shapereader import Reader import cartopy.feature as cfeat import shapefile from matplotlib.path import Path from matplotlib.patches import PathPatch import matplotlib.pyplot as plt import xarray as xr import numpy as np

第一步：读取nc文件，这里有的人用netcdf4的库，我这里用的是xarray的库。大同小异，只要能读到nc文件的变量都ok。

ncfile ='' #你的nc文件路径 lon = ncfile['XLAT'].data[0,:,:] #这里时间纬度只有1维所以第一个是0，后面是全部遍历 lat = ncfile['XLONG'].data[0,:,:] T2 = ncfile['T2'].data[0,:,:]-273.15 #这里选取了2米温度作为变量，时间为第一个

第二步：(1)掩膜，这里是根据气象家园的大佬写的代码，具体链接内容我在评论区提供：

def shp2clip (fig_cf, axes, shpfile, regionlist, pro=None): ''' 该函数用于对contourf进行掩模处理 ---关键词------------------- fig_cf : 这里指的是contourf对应的对象（不知道这样的表述是不是准确的） axes : 画图用的axes轴 shpfile : 输入shp文件 regionlist : shp文件中特定的代码，在arcmap或者是python中可以直接看到想要掩模目标的地区编号 pro : 画图对应的投影系统 ''' sf = shapefile.Reader(shpfile, encoding='gbk') vertices = [] codes = [] for shape_rec in sf.shapeRecords(): region_name = shape_rec.record[1] if region_name in regionlist: 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=axes.transData) for contour in originfig.collections: contour.set_clip_path(clip) return clip

第二步：(2)兰伯特投影刻度，有画兰伯特需求的朋友都知道，cartopy给的兰伯特投影刻度有时候比较尴尬，虽然是对的，但是为了美观才有了大佬们写的兰伯特的刻度定义函数代码：

def find_x_intersections(axes, xticks): xmin, xmax, ymin, ymax = axes.get_extent() lonmin, lonmax, latmin, latmax = axes.get_extent(ccrs.PlateCarree()) xaxis = sgeom.LineString([(xmin, ymin), (xmax, ymin)]) lon_ticks = [tick for tick in xticks if tick >= lonmin and tick = latmin and tick

【本文地址】

公司简介

联系我们