python 空间变换 python空间分布图

0.写在前面

粗略确定了一个研究范围之后，还在分析站点数据...

本文包括在中国地图上画网格和绘制站点数量空间分布图（0.5°×0.5°格子）

1.txt文件读取并提取需要的列数据

df = pd.read_table(r"D:\20140101.txt",sep='\t',header=None) 
#sep是分隔符，\t是空格，若为‘，’需要修改
#（print(df)可以查看读取的文件数据，dtype是数据格式）


lat = df[64] #df[数]是引用列的格式，我需要的纬度在第66列
#print(lat)，可以看一下有没有选错数据
del lat[0]   #我的lat的第一个数据是lat，其他的是数，所以需要把第一个删掉；之后lat的索引从1开始
lat=lat.astype(float)  #修改lat里的数据格式为float

#lon与lat类似
lon = df[66]
#print(lon)
del lon[0]
lon=lon.astype(float)

2.把站点在地图上标出来，再画上网格

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))
#ax.set_extent([70, 136, 15, 60] , crs=ccrs.PlateCarree())，这句和上面那句等价，但是我因为用这句会报错就改成上一句了

china = shpreader.Reader(r"D:\bou2_4l.dbf").geometries()
# 绘制中国国界省界九段线等等
ax.add_geometries(china, ccrs.PlateCarree(),
                  facecolor='none', edgecolor='black', zorder=1) 
ax.add_geometries(Reader(r"D:\1级河流.shp").geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor='RoyalBlue', linewidth=0.4)
plt.title('20140101站点分布图', fontsize=20,pad=20) #pad是调整标题和图之间的距离的参数

# 设置图例,shrink调整色标长度

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())

miloc = plt.MultipleLocator(0.5) #0.5意思是以0.5°为间距画线
ax.xaxis.set_minor_locator(miloc)
ax.grid(axis='x', which='both') #both可替换为major或minor，画的线不同，可自行尝试

miloc = plt.MultipleLocator(0.5)
ax.yaxis.set_minor_locator(miloc)
ax.grid(axis='y', which='both')
plt.scatter(lon, lat, c='k',s=2)
plt.savefig("D:\China2019.png", dpi=300,bbox_inches='tight')
plt.show()

作图结果：（标题没截出来）

3.（接着1）查每个格的站点数量并做图

numb=np.zeros([50,140],dtype=float,order='C')
# 创建一个全是0的二维数据，50是我区域纬度范围的2倍（因为0.5°一个格）

#假定左下为0，0
for m in range(1,1740):
    if lat[m]<55 and lat[m]>=30 and lon[m]>=70 and lon[m]<140:
        a=lat[m]
        i=int(((a-30)//0.5))
        b=lon[m]
        j=int(((b-70)//0.5))
        numb[i][j]=numb[i][j]+1
#print(numb)，可以查看现在的含数量信息的数组


plt.matshow(numb,cmap='gist_heat_r') # 生成图像
ax=plt.gca()
ax.xaxis.set_ticks_position('bottom') #把x轴刻度规定在下方
ax.invert_yaxis() #因为它作图默认(0,0)在左上，要把纵坐标颠倒一下


plt.title("2014站点数目分布图")

#plt.xticks([])

#plt.yticks([]) 需要不显示刻度的话可以加上这两句
plt.colorbar(shrink=0.8) #shrink控制colorbar长度

plt.show()

作图结果：（我的是不显示刻度的结果，需要刻度的话可以把倒数三四行加上，不过我觉得这个刻度也不是经纬度，可能没什么用）

4.把两个图叠加在一起

如果有会的兄弟姐妹可以私聊我一下吗qwq

我只能用平板的画图软件把两个图叠一起www