根据经纬度自动寻找最近的点,并返回相应属性

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

You may also like...