背景
这个是我系列插值文章的第三篇,使用机器学习插值(使用随机森林算法插值)。代码链接
代码我已经放在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/