问题描述
我正在尝试根据此处提到的过程生成数字地图，但是我不知道如何在python中实现它，尤其是对于第8步。我目前正在使用while循环进行错误测量，以确定正确的R_ref值，但要根据经验确定R_ref值的正确步长非常繁琐且困难。
有人可以建议我使用正确的函数/模块或数值方法吗？如果你们有这个的事先实现，那就太好了。
import scipy.integrate as sp
import scipy.optimize as opt
import numpy as np
import math
import random
import pandas as pd
# set source excel files
source = 'Formatted MT.xlsx'
source_1 = 'Formatted Temp.xlsx'
# set required list for reading excel files
sheet_names = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
month_rain_map_list = []
month_temp_map_list = []
days_list = [31,28.25,31,30,31]
prob_list = [0.001,0.002,0.003,0.005,0.01,0.05,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90,0.95,0.99]
#integral function of step 8
def integrand(t):
return math.exp(((-t ** 2)/2))
#take all data and store in a list of dataframes
for month in sheet_names:
month_map_rain = pd.read_excel(source,sheet_name=month,index_col=0,engine='openpyxl')
month_map_temp = pd.read_excel(source_1,engine='openpyxl')
month_rain_map_list.append(month_map_rain)
month_temp_map_list.append(month_map_temp)
#generate a dataframe for each probability in prob_list
for prob in prob_list:
dfinal = pd.DataFrame(index=month_rain_map_list[0].index,columns=month_rain_map_list[0].columns)
for lat in range(-90,91):
for lon in range(-180,181):
P_o_list = []
r_ii_list = []
r_ref =5
for i in range(len(month_rain_map_list)):
# find closest columns/index in dataframe
lon_col_rain = min(month_rain_map_list[i].columns,key=lambda x: abs(x-lon))
lat_ind_rain = min(month_rain_map_list[i].index,key=lambda x: abs(x-lat))
lon_col_temp = min(month_temp_map_list[i].columns,key=lambda x: abs(x-lon))
lat_ind_temp = min(month_temp_map_list[i].index,key=lambda x: abs(x-lat))
Mt_ii = month_rain_map_list[i].loc[lat_ind_rain].at[lon_col_rain]
N_ii = days_list[i]
t_ii = month_temp_map_list[i].loc[lat_ind_temp].at[lon_col_temp] -273.15
r_ii = 0.5874*math.exp(0.0883 * t_ii)
P_o = 100*(Mt_ii/(24*N_ii*r_ii))
if P_o > 0.70:
P_o = 0.70
r_ii = ((100/70) * ((Mt_ii)/(24*N_ii)))
r_ii_list.append(r_ii)
P_o_list.append(P_o*N_ii)
P_o_annual = sum(P_o_list)/365.25
if prob > P_o_annual:
r_ref = 0
else:
print('ATTEMPTING CONVERGENCE')
error = 1
counter = 0
while error > 0.06:
counter += 1
if error > 1:
r_ref = r_ref + (r_ref*error/2)
elif 0.60> error > 0.001:
if error > prevIoUs_error:
r_ref = r_ref + 0.001 + (r_ref*error/10)
else:
r_ref = r_ref - 0.001 - (r_ref*error/10)
else:
r_ref = r_ref - 0.01 - (r_ref*error/2)
if r_ref < 0:
r_ref = 0.01
P_ii_list = []
for i in range(len(P_o_list)):
x = (math.log(r_ref) + 0.7938 - math.log(r_ii_list[i]))/1.26
q_x = (1/(math.sqrt(2*math.pi))) * sp.quad(integrand,x,np.inf)[0]
P_ii = P_o_list[i] * q_x
P_ii_list.append(P_ii)
P_o_final = sum(P_ii_list)/365.25
prevIoUs_error = error
error = abs(P_o_final/prob - 1)
if r_ref < 0.005 and error > 0.05:
r_ref = 0
break
if r_ref < 1 and P_o_final < 0.00001:
r_ref = 0
break
if counter > 10000:
print('max limit hit')
print(f'r_ref:{r_ref},\n P_o_final = {P_o_final},\n prob: {prob},\n error: {error}')
print(f' P_ii_list : {P_ii_list}')
raise Exception # generally the code breaks here as it cannot hit convergence
print('CONVERGED')
print(f'lat: {lat},lon: {lon},Mt_ii = {Mt_ii},r_ii = {r_ii},t_ii = {t_ii},P_o = {P_o_list}')
print(f'r_ref: {r_ref}')
dfinal.loc[lat_ind_rain].at[lon_col_rain] = r_ref
print(dfinal)
raise Exception
print(dfinal)
print(month_map_list[0])
我的代码的最终目标是为prob_list中定义的每种超出概率生成降雨率图。国际电联只提供0.01％的地图。
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）