python降水 | 自定义挑选部分站点判断代表性

0.写在前面

任务二：找到一种选取5°×5°格子内部分站点的方法，挑选站点，并判断站点是否有代表性

难点：如何选站点

解决办法：将每个方格按九宫格划分，选位于四角和中心的站点。

1.完整代码

from importlib.resources import path
import os
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from cartopy.io.shapereader import Reader
from cartopy.io import shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.crs as ccrs
import time
import xarray as xr
time_start = time.time()
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

def dave(df):
    lat = df[64]
    del lat[0]
    lat=lat.astype(float)
    lon = df[66]
    del lon[0]
    lon=lon.astype(float)
    pre=df[72]
    del pre[0]
    pre=pre.astype(float) 
    daypre=np.zeros([5,14],dtype=float,order='C')
    daysta=np.zeros([5,14],dtype=float,order='C') 
    daypre9=np.zeros([5,14],dtype=float,order='C') 
    daysta9=np.zeros([5,14],dtype=float,order='C')
    for m in range(1,(len(lat)+1)): 
        if lat[m]<55 and lat[m]>=30 and lon[m]>=70 and lon[m]<140:
            a=lat[m]
            i=4-int(((a-30)//5))
            e=int((a-30)/1.666667)
            b=lon[m]
            j=int(((b-70)//5))
            r=int((b-70)/1.666667)
            if pre[m]<500:
                daysta[i][j]=daysta[i][j]+1
                daypre[i][j]=daypre[i][j]+pre[m]
            str1={0,2,3,5,6,8,9,11,12,14}
            str2={0,3,6,9,12,15,18,21,24,27,30,33,36,39,2,5,8,11,14,17,20,23,26,29,32,35,38,41}
            q=(((e in str1) and (r in str2))or((not(e in str1)) and (not(r in str2))))
            if pre[m]<500 and q:  #判断是否在要选的格子里且没有缺测
                daysta9[i][j]=daysta9[i][j]+1
                daypre9[i][j]=daypre9[i][j]+pre[m]
    dayave=np.zeros([5,14],dtype=float,order='C') 
    dayave9=np.zeros([5,14],dtype=float,order='C')
    for i in range(0,5):
        for j in range(0,14):
            if daysta[i][j]!=0:
                dayave[i][j]=daypre[i][j]/daysta[i][j]
            if daysta9[i][j]!=0:
                dayave9[i][j]=daypre9[i][j]/daysta9[i][j]
    return dayave,dayave9



def draw(ave):
    fig = plt.figure(1, figsize=[16, 9])
    proj = ccrs.PlateCarree()
    ax = plt.subplot(1, 1, 1, projection=proj)
    extent = [70, 140, 30, 55]
    ax.set(xlim=(70,140),ylim=(30,55))
    china = shpreader.Reader(r"D:\bou2_4l\bou2_4l.dbf").geometries()
    plt.title('2014降水总量相对误差', fontsize=20,pad=20)
    ax.add_geometries(china, ccrs.PlateCarree(),
                  facecolor='none', edgecolor='black', zorder=2) 
    ax.add_geometries(Reader(r"D:\1~5级水系矢量图\river\1级河流.shp").geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='RoyalBlue', linewidth=0.4)             
    c=np.ones([5,14],dtype=float,order='C')
    alpha=np.where(ave!=0,c*0.7,0)
    plt.imshow(ave,cmap='rainbow',extent=(70,140,30,55),alpha=alpha,zorder=3)
    cb=plt.colorbar(shrink=0.5)
    ax.set_xticks(np.arange(extent[0], extent[1]+1, 10), crs=proj)
    ax.set_yticks(np.arange(extent[-2], extent[-1]+1, 5), crs=proj)
    ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
    ax.yaxis.set_major_formatter(LatitudeFormatter())
    plt.tick_params(labelsize=12)
    plt.show()
 

path = r"D:\2014-2020\2014"
files= os.listdir(path)
k=1
yearsum=np.zeros([5,14],dtype=float,order='C')
ninesum=np.zeros([5,14],dtype=float,order='C')
for filename in files:
    df = pd.read_table(os.path.join(path,filename),sep='\t',header=None)
    day,sele=dave(df)
    yearsum=yearsum+day
    ninesum=ninesum+sele

alpha=np.zeros([5,14],dtype=float,order='C')  #alpha 是相对误差
for i in range(0,5):
        for j in range(0,14):
            if yearsum[i][j]!=0:# and ninesum[i][j]!=0:
                alpha[i][j]=(ninesum[i][j]-yearsum[i][j])/yearsum[i][j]

draw(alpha)

#筛选不到站点的地方值是-1
#0：除了没有站点的地方，还有站点全部被选到的格子

time_end = time.time()  # 记录结束时间
time_sum = time_end - time_start  # 计算的时间差为程序的执行时间，单位为秒/s
print('运行时间:',time_sum)

2. 成图：