python 气象 | 白化、南海子图

最新推荐文章于 2024-08-09 15:34:32 发布

记事本本本本

最新推荐文章于 2024-08-09 15:34:32 发布

阅读量540

点赞数 4

分类专栏：气象python实习报告文章标签： python

本文链接：https://blog.csdn.net/qq_51633944/article/details/124956236

版权

短期气候预测实习报告 python 中国雨型

以前的东西备份在这里

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.mpl.ticker as cticker
from cartopy.util import add_cyclic_point
from cartopy.io.shapereader import Reader
import os
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
# from wrf import smooth2d
from scipy import stats
from scipy.stats import ttest_1samp
from xgrads import open_CtlDataset
#maskout----------------------------
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def shp2clip(originfig,ax,shpfile,region):
    sf = shapefile.Reader(shpfile)
    for shape_rec in sf.shapeRecords():
        if shape_rec.record[3] == region:  ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。
            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 += [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

SHP = r"D:\【短期气候预测】\china_shp"
def make_map(ax,box):
    # 加国界 
    ax.add_geometries(Reader(os.path.join(SHP, 'country1.shp')).geometries(),
                      ccrs.PlateCarree(),facecolor='none',edgecolor='k', linewidth=0.5)
       #标注坐标轴
    gl=ax.gridlines(draw_labels=True,linestyle=":",linewidth=0.1 ,x_inline=False, y_inline=False,color='k')
    gl.top_labels=False; gl.right_labels=False #关闭上部和右部经纬标签 
    gl.rotate_labels = None  # 禁止gridlines旋转正向
    gl.xformatter = LONGITUDE_FORMATTER  #使横坐标转化为经纬度格式            
    gl.yformatter = LATITUDE_FORMATTER                                        
    gl.xlocator=mticker.FixedLocator(np.arange(box[0],box[1],5))                                   
    gl.ylocator=mticker.FixedLocator(np.arange(box[2],box[3],5)) 
    gl.xlabel_style={'size':10}; gl.ylabel_style={'size':10}#修改经纬度字体大小                                 
    ax.spines['geo'].set_linewidth(0.5)#调节边框粗细
    ax.set_extent([box[0],box[1],box[2],box[3]],crs=ccrs.PlateCarree())
    #添加网格线
    ax.gridlines(linestyle='--',alpha=0.4)   
    return ax
def plot_ax(i,ax,lon,lat,var,proj,ccmap,lable=False):
    ax = make_map(ax,box)
    # if ifsmooth:
    #     var = smooth2d(var, 3, cenweight=4) #  平滑数据  3：平滑次数 cenweight=4 平滑权重
    # contourf
    cf = ax.contourf(lon,lat,var,cmap = cc