Python | 中国站点数量分布图（0.25°小方格）

0.写在前面

之前写过一版画图代码，但是matshow函数不太好用，这次改为imshow函数绘制

1.数据准备

计算每个小方格中的站点数量

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import time
time_start = time.time()

#计算30-55N，70-140E的2015、2017年站点数目
#0.25°的格子


def stas(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) #999999不算降水，0算
    numb=np.zeros([120,280],dtype=float,order='C') 
    for m in range(1,len(lat)+1): #遍历一天内每个站点
        if lat[m]<55 and lat[m]>=25 and lon[m]>=70 and lon[m]<140:
            a=lat[m]
            i=119-int(((a*100-2500)//25))
            b=lon[m]
            j=int(((b*100-7000)//25))
            numb[i][j]=numb[i][j]+1
    return numb

sta1=np.zeros([120,280],dtype=float,order='C')
sta2=np.zeros([120,280],dtype=float,order='C')
plus=np.zeros([120,280],dtype=float,order='C')
df1 = pd.read_table(r"E:\2014-2020\2015\SURF_CHN_MUL_DAY_20150101.txt",sep='\t',usecols=(64,66,72),low_memory=False,header=None)
sta1=stas(df1)
df2 = pd.read_table(r"E:\2014-2020\2017\SURF_CHN_MUL_DAY_20171231.txt",sep='\t',usecols=(64,66,72),low_memory=False,header=None)
sta2=stas(df2)
plus=sta2-sta1

np.savetxt(r"D:\20150101站点数目.txt", sta1, delimiter ="\t")
np.savetxt(r"D:\20171231站点数目.txt", sta2, delimiter ="\t")
time_end = time.time()  # 记录结束时间
time_sum = time_end - time_start  # 计算的时间差为程序的执行时间，单位为秒/s
print('运行时间:',time_sum)

2.画图

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
time_start = time.time()
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
sta1 = np.loadtxt(r"D:\20171231站点数目.txt")
sta= np.loadtxt(r"D:\20150101站点数目.txt")

fig = plt.figure(figsize=(10,10),dpi=300)#定义图片像素
proj = ccrs.PlateCarree()
extent = [70, 140, 25, 55]
china = shpreader.Reader(r"D:\bou2_4l\bou2_4l.dbf").geometries()
c=np.ones([120,280],dtype=float,order='C')


ax = plt.subplot(2, 1, 1, projection=proj)
ax.set(xlim=(70,140),ylim=(25,55))
plt.title('(a)2015年站点数目',loc='left',fontsize=20)# 绘制标签
ax.add_geometries(china, ccrs.PlateCarree(),
                  facecolor='none', edgecolor='black', zorder=2,linewidth=0.5) 
ax.add_geometries(Reader(r"D:\river\1级河流.shp").geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='RoyalBlue', linewidth=0.4)             
alpha=np.where(sta!=0,c*0.8,0)
plt.imshow(sta,cmap='rainbow',extent=(70,140,25,55),alpha=alpha,zorder=3,vmin=0,vmax=40)
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=15)



#____________________________________________________________________________________________________

ax1 = plt.subplot(2, 1, 2, projection=proj)
ax1.set(xlim=(70,140),ylim=(25,55))
plt.title('(b)2017年站点数目',loc='left',fontsize=20)# 绘制标签
china = shpreader.Reader(r"D:\bou2_4l\bou2_4l.dbf").geometries()
ax1.add_geometries(china, ccrs.PlateCarree(),
                  facecolor='none', edgecolor='black', zorder=2,linewidth=0.5) 
ax1.add_geometries(Reader(r"D:\river\1级河流.shp").geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='RoyalBlue', linewidth=0.4)             
alpha=np.where(sta1!=0,c*0.8,0)
plt.imshow(sta1,cmap='rainbow',extent=(70,140,25,55),alpha=alpha,zorder=3,vmin=0,vmax=40)
ax1.set_xticks(np.arange(extent[0], extent[1]+1, 10), crs=proj)
ax1.set_yticks(np.arange(extent[-2], extent[-1]+1, 5), crs=proj)
ax1.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax1.yaxis.set_major_formatter(LatitudeFormatter())
plt.tick_params(labelsize=15)


l = 0.16
b = 0.06
w = 0.7
h = 0.02

# 对应 l,b,w,h；设置colorbar位置；
rect = [l, b, w, h]

cbar_ax = fig.add_axes(rect)
cb=plt.colorbar(shrink=1,extend='both',cax=cbar_ax,orientation='horizontal')
cb.ax.tick_params(labelsize=15)
cb.set_label('个',fontsize=15)
#——————————————————————————————————————————————————————————————————————————————————————————————————————————————



plt.savefig(r"D:\站点数目.png", dpi=300)
plt.show()

3.出图

 

ps.不要在意数据是否符合实际，仅作为例子展示