背景
这个是我系列插值文章的第三篇，使用机器学习插值（使用随机森林算法插值）。代码链接
代码我已经放在Github上面了，免费分享使用，/yuanzhoulvpi/tiny_python/tree/main/python_GIS。
介绍
本文是python与GIS数据处理系列中的插值部分————使用机器学习算法插值（随机森林算法插值）。
我这里的方法并不是最简单的方法，并不是一行代码就能从头到尾实现这个插值功能的，我这里的目的是：使用python完善的数学工具链，从想法到代码实现，一套完整的数据处理思路。
过程
1. 准备数据
历史数据
我们插值，需要利用特定的历史数据，才能做插值。不然依靠的数据基础是什么？
import numpy as npimport pandas as pdimport matplotlib.pyplot as plt# part 1name_list = ['id', 'latitude', 'longitude', 'altitude', 'year', 'month', 'day', 'avg_temp', 'max_temp','min_temp', 'col1', 'col2', 'col3']raw_data = pd.read_table("./数据集/插值数据/SURF_CLI_CHN_MUL_DAY-TEM-12001-05.TXT",header=None, sep='\t')raw_data = raw_data[0].apply(lambda x: pd.Series([float(i.strip()) for i in x.split(' ') if i != ''], index=name_list))clean_data = raw_data.loc[(raw_data['year'] == ) & (raw_data['month'] == 5), :].groupby(['latitude', 'longitude']).agg(avg_temp=('avg_temp', 'mean')).reset_index()# part 2for index in ['latitude', 'longitude']:clean_data[index] = clean_data[index] * 0.01clean_data = clean_data.loc[(clean_data['longitude'] > 20) & (clean_data['latitude'] > 10), :]clean_data
上面的一番操作，具体做的内容是在一个时间段下，对数据做处理，清洗，大概得出一个数据框，这个数据库主要包括：经纬度和变量，这里的变量就是avg_temp。
确定地理边界
插值，也不是盲目的插值，要选定好特定的区域，既然是安徽人，那我就直接选择安徽省作为边界。
from getchinamap import getchinamapchina_engine = getchinamap.DownloadChmap()prov_gpd = china_engine.download_province(province_name='安徽省', target='边界')# prov_gpd = china_engine.download_country(target='边界')prov_gpd
上面那个包是我写的，使用高德数据和简单的代码，就可以获得完整的、准确的中国各级地图数据。
2. 数据处理
数据获得后，我们要对数据做过大概的估计，对点做一些大概的筛选（筛选在安徽省的点），还要对点做一些预处理之类的。
prov_gpd_valid = prov_gpd.copy()prov_gpd_valid['geometry'] = prov_gpd_valid.buffer(0)from shapely.geometry import Pointdef detect_pic(x):return prov_gpd_valid.contains(Point(x['longitude'], x['latitude']))[0]def detect_pic_inrect(x):bounds_ = prov_gpd_valid.bounds.iloc[0, :]minx, miny, maxx, maxy = bounds_.minx, bounds_.miny, bounds_.maxx, bounds_.maxyreturn (minx <= x['longitude']) & (x['longitude'] <= maxx) & (x['latitude'] >= miny) & (x['latitude'] <= maxy)# clean_data['in_geo'] = clean_data.apply(lambda x: detect_pic(x), axis=1)clean_data['in_box'] = clean_data.apply(lambda x: detect_pic_inrect(x), axis=1)prov_pointer_df = clean_data.loc[clean_data['in_box']].reset_index(drop=True)prov_pointer_df
可以看出来，大概有37个点在安徽境内。
然后再画图看看。
fig, ax = plt.subplots(figsize=(4, 6))ax.scatter(prov_pointer_df['longitude'], prov_pointer_df['latitude'])# ax.scatter(clean_data['longitude'], clean_data['latitude'])prov_gpd_valid.boundary.plot(ax=ax)ax.set_xlim(prov_gpd_valid.bounds.minx[0], prov_gpd_valid.bounds.maxx[0])ax.set_ylim(prov_gpd_valid.bounds.miny[0], prov_gpd_valid.bounds.maxy[0])
然后再看看安徽省的上下左右的具体数值。
3. 插值处理
# part 1# 将安徽省的边界提取出来，用来生成网格bounds_ = prov_gpd_valid.bounds.iloc[0, :]minx, miny, maxx, maxy = bounds_.minx, bounds_.miny, bounds_.maxx, bounds_.maxylongitude_x = np.linspace(start=minx, stop=maxx, num=100)latitude_y = np.linspace(start=miny, stop=maxy, num=200)# part 2# 使用历史数据做插值的函数拟合。这里选择的是quintic方法from sklearn.ensemble import RandomForestRegressorx_train = clean_data[['latitude', 'longitude']]y_train = clean_data[['avg_temp']]# 这里需要调整模型的参数rf_model = RandomForestRegressor() #n_estimators=2000,min_samples_split=20rf_model.fit(x_train, y_train)# part 3# 使用拟合的函数做预测predict_Latitude, predict_longitude = np.meshgrid(latitude_y, longitude_x)predict_df = pd.DataFrame({'latitude':predict_Latitude.flatten(),'longitude':predict_longitude.flatten()})predict_value = rf_model.predict(predict_df)predict_matrix = predict_value.reshape(predict_Latitude.shape)predict_matrix.shape
关于上面的插值部分，要详细介绍一下：
part1 部分就是为了将安徽省的边界（四个角）都提取出来，然后把横向的、纵向的都生成稠密的点，经纬度都稠密的话，不就是变成网格了么。这个时候，你需要仔细观察longitude_x和latitude_y的shape。如果觉得100， 200 太大了，就改成10， 20，看看。part2 部分就是将历史历史数据放在了函数中，因为经纬度是二维的，肯定是使用2d。传递的参数，分别是x， y，z。这里的x，y，z都是上面看到的面板数据（就是各个点的数据),不是什么网格数据。在part2中，这里的随机森林的参数是没有调整的，如果想调整，可以看看官网了解一下：https://scikit-/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html在part3 部分，返回的是网格数据，也就是一个二维的预测数据。
4. 可视化
# part 1import matplotlib as mplfrom matplotlib import cmfig, ax = plt.subplots(figsize=(10, 6), dpi=150)colors = ["#33A02C", "#B2DF8A", "#FDBF6F", "#1F78B4", "#999999", "#E31A1C"]# ax.scatter(clean_data['longitude'], clean_data['latitude'], c='red')prov_gpd_valid.boundary.plot(ax=ax, color='white')# ax_im_bar = ax.contourf(grid_x, grid_y, predict_cubic,cmap=cm.coolwarm) #mpl.colors.LinearSegmentedColormap.from_list("mypalette", colors, N=1000)# part 2# 将刚才生成的预测数值转换成图像并且画出来ax_im_bar = ax.imshow(predict_matrix, origin='lower',extent=(minx, maxx,miny, maxy),cmap=mpl.colors.LinearSegmentedColormap.from_list("mypalette", colors, N=1000))# ax.contour(grid_x, grid_y, predict_cubic)ax.scatter(prov_pointer_df['longitude'], prov_pointer_df['latitude'], c='black', s=6)# part 3# 画出历史数据的点出来for index in range(prov_pointer_df.shape[0]):ax.text(prov_pointer_df.iloc[index]['longitude'], prov_pointer_df.iloc[index]['latitude'],np.around(prov_pointer_df.iloc[index]['avg_temp'], 2), c='black')## ax.set_xlim(prov_gpd_valid.bounds.minx[0], prov_gpd_valid.bounds.maxx[0])# ax.set_ylim(prov_gpd_valid.bounds.miny[0], prov_gpd_valid.bounds.maxy[0])# part 4 # 显示color barfig.colorbar(ax_im_bar, orientation='vertical')plt.tight_layout()plt.savefig("结果/result022701.png")
参考
查看我的历史仓库代码：/yuanzhoulvpi/tiny_python/tree/main/python_GIS
了解更多的随机森林算法细节，可以查看这个链接：https://scikit-/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html
查看下载中国地图数据包：/project/getchinamap/
