根据经纬度自动寻找最近的点,并返回相应属性
import pandas as pd
import numpy as np
from scipy.spatial.distance import cdist
import scipy.stats as stats
from pylab import *
STA_tmin = 'F:/AAA/STA-T-TMAX-TMIN-CSV/STA-TMEAN.csv'
PG_tmin = 'F:/AAA/GRID-TEM-CSV/IF-MEAN/PGTMEAN.csv'
sta_tmin = pd.read_csv(STA_tmin)
pg_tmin = pd.read_csv(PG_tmin)
# 两组数据
def closest_point(point, points): # 根据哪列,在哪里找
""" Find closest point from a list of points. """
return points[cdist([point], points).argmin()] # 前者必须是个array
def match_value(df, col1, x, col2): # 哪里找,根据什么找,和谁匹配,返回什么值。
""" Match value x from col1 row to value in col2. """
return df[df[col1] == x][col2].values[0]
pg_tmin['point'] = [(x, y) for x, y in zip(pg_tmin['lat'], pg_tmin['lon'])]
sta_tmin['point'] = [(x, y) for x, y in zip(sta_tmin['lat'], sta_tmin['lon'])]
pg_tmin = (pg_tmin.groupby('point').mean()).reset_index()
sta_tmin = (sta_tmin.groupby('point').mean()).reset_index()
sta_tmin['closest'] = [closest_point(x, list(pg_tmin['point'])) for x in sta_tmin['point']] # 遍历point 找出离得最近得点,双双验证
pg_tmin['closest'] = [closest_point(x, list(sta_tmin['point'])) for x in pg_tmin['point']] # 遍历point # 找出离得最近得点,双双验证
p = (sta_tmin['closest']).values # 和实测靠近的pg点 127个
p2 = (pg_tmin['closest']).values # 和pg靠近得实测点 134个
pg_tmin = pg_tmin[pg_tmin['point'].isin(p)] # 筛选出pg和实测靠近得点,closet自动丢掉了某些很近的实测点,再用筛出的实测中再从实测筛选一遍。
pg_tmin.drop_duplicates('closest','first',inplace = True) # 去重
pg_tmin = pg_tmin[(pg_tmin['year']>=1970)&(pg_tmin['year']<=1989)] # 拿出pg和实测一样年份的数据
sta_tmin = sta_tmin[sta_tmin['point'].isin((pg_tmin['closest']).values)] # 再用筛出的实测中再从实测筛选一遍。
pg_tmin['ele-sta'] = [match_value(sta_tmin, 'point', x, 'ele') for x in pg_tmin['closest']] # 返回到pg对应得实测站点海拔
pg_tmin.sort_values(by ='ele-sta',inplace=True) # 按照实测站点海拔排序
sta_tmin.sort_values(by = 'ele',inplace=True) # 按照实测站点海拔排序
pg_tmin = pg_tmin.reset_index()
sta_tmin = sta_tmin.reset_index()
ele_sta = pg_tmin['ele-sta'] # 实测海拔
ele_pg = pg_tmin['ele-grid']
pg_tmin['sub']= pg_tmin['ele-grid']- pg_tmin['ele-sta'] # 海拔差
months = ['one', 'two', 'three', 'four', 'five', 'six', 'seven', 'eight', 'nine','ten', 'eleven', 'december']
for month in months:
pg_tmin[month]=pg_tmin[month]+pg_tmin['sub']*0.7/100 # 移动到同一点
plt.figure()
plt.scatter(ele_sta,sta_tmin['three'])
x = ele_sta
y = (sta_tmin['three']).values
z = np.polyfit(x, y, 1)
p = np.poly1d(z)
r, p_value=stats.pearsonr(x,y)
plt.plot(x, p(x), '--k', dashes=(2, 14))
plt.scatter(ele_sta,pg_tmin['three'])
plt.show()
import matplotlib.pyplot as plt # 画每个月的图
for month in months:
plt.figure()
minorticks_on()
tick_params(which='major', width=2)
ax = plt.gca()
ax.spines['top'].set_visible(False) # 去掉上边框
ax.spines['right'].set_visible(False)
x = (pg_tmin[month]).values
y = (sta_tmin[month]).values
plt.scatter(pg_tmin[month],sta_tmin[month],color='black',s = 33)
plt.xticks(fontsize = 18,)
plt.yticks(fontsize = 18,)
z = np.polyfit(x, y, 1)
p = np.poly1d(z)
r, p_value=stats.pearsonr(x,y)
plt.plot(x, p(x), '--k', dashes=(2, 14))
junk = str((month,round(z[0],2),round(z[1],2),round(r*r,2)))
print(month,round(z[0],2),round(z[1],2),round(r*r,2),'p=',p_value)
# plt.savefig('F:/PPPPSSSSS/TMIN-LINE/'+month+'-tmin'+'.tif',dpi=600,bbox_inches='tight',pad_inches = 0)
plt.show()
# with open('F:/PPPPSSSSS/TMAX-FORMULA.txt','at') as f:
# print(*junk, sep=' ', file=f)
转载自:https://blog.csdn.net/hengcall/article/details/82807946