python与GIS数据处理——随机森林算法插值

背景

1. 这个是我系列插值文章的第三篇，使用机器学习插值（使用随机森林算法插值）。

代码链接

代码我已经放在Github上面了，免费分享使用，https://github.com/yuanzhoulvpi2017/tiny_python/tree/main/python_GIS。

介绍

1. 本文是python与GIS数据处理系列中的插值部分————使用机器学习算法插值（随机森林算法插值）。

2. 我这里的方法并不是最简单的方法，并不是一行代码就能从头到尾实现这个插值功能的，我这里的目的是：使用python完善的数学工具链，从想法到代码实现，一套完整的数据处理思路。

过程

1. 准备数据

历史数据

我们插值，需要利用特定的历史数据，才能做插值。不然依靠的数据基础是什么？

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# part 1
name_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-201905.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'] == 2019) & (raw_data['month'] == 5), :].groupby(
    ['latitude', 'longitude']).agg(
    avg_temp=('avg_temp', 'mean')
).reset_index()

# part 2
for index in ['latitude', 'longitude']:
    clean_data[index] = clean_data[index] * 0.01
clean_data = clean_data.loc[(clean_data['longitude'] > 20) & (clean_data['latitude'] > 10), :]
clean_data

上面的一番操作，具体做的内容是在一个时间段下，对数据做处理，清洗，大概得出一个数据框，这个数据库主要包括：经纬度和变量，这里的变量就是avg_temp。

确定地理边界

插值，也不是盲目的插值，要选定好特定的区域，既然是安徽人，那我就直接选择安徽省作为边界。

from getchinamap import getchinamap

china_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 Point


def 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_.maxy
    return (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_.maxy

longitude_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 RandomForestRegressor
x_train = clean_data[['latitude', 'longitude']]
y_train = clean_data[['avg_temp']]

# 这里需要调整模型的参数
rf_model = RandomForestRegressor() #n_estimators=2000,min_samples_split=20
rf_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

关于上面的插值部分，要详细介绍一下：

1. part1 部分就是为了将安徽省的边界（四个角）都提取出来，然后把横向的、纵向的都生成稠密的点，经纬度都稠密的话，不就是变成网格了么。这个时候，你需要仔细观察longitude_x和latitude_y的shape。如果觉得100， 200 太大了，就改成10， 20，看看。
2. part2 部分就是将历史历史数据放在了函数中，因为经纬度是二维的，肯定是使用2d。传递的参数，分别是x， y，z。这里的x，y，z都是上面看到的面板数据（就是各个点的数据),不是什么网格数据。
3. 在part2中，这里的随机森林的参数是没有调整的，如果想调整，可以看看官网了解一下：https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html
4. 在part3 部分，返回的是网格数据，也就是一个二维的预测数据。

4. 可视化

# part 1
import matplotlib as mpl
from matplotlib import cm


fig, 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 bar
fig.colorbar(ax_im_bar, orientation='vertical')
plt.tight_layout()
plt.savefig("结果/result022701.png")

参考

1. 查看我的历史仓库代码：https://github.com/yuanzhoulvpi2017/tiny_python/tree/main/python_GIS

2. 了解更多的随机森林算法细节，可以查看这个链接：https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html

3. 查看下载中国地图数据包：https://pypi.org/project/getchinamap/