考虑海拔的IDW的插值_idw结合dem插值

import pandas as pd
from geopy.distance import vincenty
import numpy as np
loworhigh = 'high'  # 输入需要哪段 ,low 代表第一段，high代表第二段。
year = 1951     # 是否要选择年份
dem = pd.read_table('F:/AAA/PG-T-LAPSE-COR/lldem_syr',header=None)      # 流域模拟点高程文件
sta_df = pd.read_csv('F:/AAA/P-LAPSE/STA-SYR-PRECP-data-interpolate.csv') # 降水文件
sta_df = sta_df[sta_df['year']==year]
sta_df['point'] = [(x, y) for x, y in zip(sta_df['lat'], sta_df['lon'])]
sta_df = sta_df.groupby('point').mean().reset_index()
 
dem.columns = ['lat','lon','ele']
dem['point'] = [(x, y) for x, y in zip(dem['lat'], dem['lon'])]
 
lapse_low = [0.057,0.057,0.078,0.08,0.058,0.021,0.019,0.013,0.008,0.063,0.067,0.074]         # 插值梯度1
lapse_high = [-0.012,-0.014,-0.022,-0.023,-0.01,0.021,0.019,0.013,0.008,-0.008,-0.01,-0.01]  # 插值梯度2
 
dem_low = dem[dem['ele']<=1256]   # 拿出小于1300的模拟格子，用第一段插值。
dem_high = dem[dem['ele']>1256]
 
sta_df_low = sta_df[sta_df['ele']<=1256] # 拿出小于1256的实测点
sta_df_high = sta_df[sta_df['ele']>1256]  # 拿出大于1256的实测点
 
 
months = ['one','two','three','four','five','six','seven','eight','nine','ten','eleven','december']
months_latlon = ['lat','lon','one','two','three','four','five','six','seven','eight','nine','ten','eleven','december']
a = []
b = []
junk = pd.DataFrame()
result_annual = pd.DataFrame(np.zeros((2,14)),columns=months_latlon)
if loworhigh=='high':
    sta_df_run = sta_df_high
    dem_run = dem_high
    dem_run = dem_run.reset_index()
    lapse_run = lapse_high
else:
    sta_df_run =sta_df_low
    dem_run = dem_low
    dem_run = dem_run.reset_index()
    lapse_run = lapse_low
 
for point in dem_run['point'].values:
    for point_sta in sta_df_run['point'].values:
        dist = vincenty(point, point_sta).km
        print(point)
        print(point_sta)
        a.append(dist)
        b.append(point_sta)
    lat = point[0]
    lon = point[1]
    grid_ele = dem_run['ele'][(dem_run['lat']==lat)&(dem_run['lon']==lon)]  # 预测点所在海拔
    junk['distance'] = a
    junk['sta_point'] = b
    junk.sort_values(by = 'distance',inplace=True)
    nearst_point = junk['sta_point'][0:3].values  # 选取附近三个实测点，根据实际需要调整
    sta_out_point = sta_df_run[sta_df_run['point'].isin(nearst_point)].reset_index()  # 拿出实测点，准备计算此网格
    sta_out_point['ele_grid'] = int(grid_ele)
    a = []
    b = []
    sta_out_point['sub'] = (sta_out_point['ele_grid'] - sta_out_point['ele'])
    sta_out_point['out1'] = sta_out_point['one']+sta_out_point['sub']*lapse_run[0]
    sta_out_point['out2'] = sta_out_point['two']+sta_out_point['sub']*lapse_run[1]
    sta_out_point['out3'] = sta_out_point['three']+sta_out_point['sub']*lapse_run[2]
    sta_out_point['out4'] = sta_out_point['four']+sta_out_point['sub']*lapse_run[3]
    sta_out_point['out5'] = sta_out_point['five']+sta_out_point['sub']*lapse_run[4]
    sta_out_point['out6'] = sta_out_point['six']+sta_out_point['sub']*lapse_run[5]
    sta_out_point['out7'] = sta_out_point['seven']+sta_out_point['sub']*lapse_run[6]
    sta_out_point['out8'] = sta_out_point['eight']+sta_out_point['sub']*lapse_run[7]
    sta_out_point['out9'] = sta_out_point['nine']+sta_out_point['sub']*lapse_run[8]
    sta_out_point['out10'] = sta_out_point['ten']+sta_out_point['sub']*lapse_run[9]
    sta_out_point['out11'] = sta_out_point['eleven']+sta_out_point['sub']*lapse_run[10]
    sta_out_point['out12'] =sta_out_point['december'] +sta_out_point['sub']*lapse_run[11]
    sta_out_point['lat_grid'] = point[0]
    sta_out_point['lon_grid'] = point[1]
    sta_out_point = sta_out_point[['lat_grid','lon_grid','out1','out2','out3','out4','out5','out6','out7','out8','out9','out10','out11','out12']]
    sta_out_point[sta_out_point <= 0] = 1
    sta_out_point.columns = ['lat','lon','one','two','three','four','five','six','seven','eight','nine','ten','eleven','december'] # 输出习惯的名字
    a = []
    b = []
 
    f = pd.concat([result_annual,sta_out_point],axis=0,ignore_index=True)
    result_annual = f
 
f = f.groupby(['lat','lon']).mean().reset_index(0)
f.to_csv('F:/AAA/P-LAPSE/inter-monthly2.csv') # 输出文件