griddata三维空间插值
从这一篇文章,你将要学到
- 如何利用griddata进行三维空间插值;
- 及其适用范围和进阶的逐步插值
背景
最近在做一个项目,要为上海市13000+个普通住宅楼盘算基本价格,俗称基价,可以从第三方来的案例数据只能覆盖大约3000个楼盘,余下的10000楼盘难为无米之炊,联想到地形图的思想,把上海市所有楼盘的基价看成海拔,楼盘的经纬度就是位置所在,然后会在三维空间形成一个连续平滑的三维曲面,这里利用scipy的interpolate类里面的griddata函数小试牛刀。
数据
从原数据我们看到需要插值的thismonthprice有大量空缺,如何利用地理位置进行插值呢?基本思路如下
- 将数据分成两部分,一部分是thismonthprice有价格,一部分是thismonthprice为空的;
- 画出有价格和没价格的楼盘散点图,方便直观感受;
- 利用thismonthprice有价格的进行三维曲面建模训练;
- 利用训练好的模型对thismonthprice为空的进行模拟插值。
完整代码
import numpy as np #导入数值计算模块
import pandas as pd #导入数据分析模块
from scipy import stats
import matplotlib.pyplot as plt #导入绘图模块
from mpl_toolkits.mplot3d import Axes3D #导入3d绘图模块
from matplotlib import cm #颜色调整用
from scipy.interpolate import griddata #导入插值类def scatterPlot(data): #绘制绘图函数plt.figure(figsize = (9, 7)) #新建画布caseData = data.dropna(subset = ["thismonthprice"]) #去掉空值,保留本月有价格的数据留作插值依据lon = caseData["lon"] #经度lat = caseData["lat"] #纬度interpolateData = data[data["thismonthprice"].isnull()] #待插值的数据x = interpolateData["lon"] #待插值的楼盘经度y = interpolateData["lat"] #待插值的楼盘纬度plt.scatter(lon, lat, s = 5, color = "r", label = "no need interpolate") #不需要插值的散点图plt.scatter(x, y, s = 5, color = "b", alpha = 0.2, label = "need interpolate") #需要插值的散点图plt.legend()plt.show()def gridInterpolate(data): #定义插值函数,data需要有newdiskid, lon,lat, 上个月价格,本月价格(有空缺)casenum = [] #用来存放有价格的楼盘数interpolatenum = [] #用来存放待插值的样本数caseData = data.dropna(subset = ["thismonthprice"]) #去掉空值,保留本月有价格的数据留作插值依据casenum.append(len(caseData)) #追加有价格的样本数casediskid = caseData["newdiskid"] #楼盘idlon = caseData["lon"] #经度lat = caseData["lat"] #纬度z = caseData["thismonthprice"] #价格interpolateData = data[data["thismonthprice"].isnull()] #待插值的数据interpolatenum.append(len(interpolateData)) #追加待插值样本数interpolatediskid = interpolateData["newdiskid"] #待插值的楼盘idx = interpolateData["lon"] #待插值的楼盘经度y = interpolateData["lat"] #待插值的楼盘纬度xx, yy = np.meshgrid(x, y) #平面网格点待插值的位置zz = griddata((lon, lat), z, (xx, yy), method='nearest', fill_value = z.median(), rescale = True) #stats.mode(z)[0][0]在没有价格的网格点插值Z = pd.DataFrame(zz) #数据框化interpolateprice = [] #用来存放插值for i in range(Z.shape[0]):interpolateprice.append(Z[i][i]) #取对角线上的元素label1 = pd.Series(np.zeros(len(casediskid))) #有价格的标签casedf = pd.concat([casediskid, z], axis = 1).reset_index()[["newdiskid", "thismonthprice"]] #索引重置casedf = pd.concat([casedf, label1], axis = 1)casedf.columns = ["newdiskid", "z", "label"]label2 = pd.Series(np.ones(len(interpolatediskid))) #插值标签interpolatedf = pd.DataFrame({"interpolatediskid":interpolatediskid, "interpolateprice": interpolateprice}).reset_index()[["interpolatediskid", "interpolateprice"]]interpolatedf = pd.concat([interpolatedf, label2], axis = 1)interpolatedf.columns = ["newdiskid", "z", "label"] result = pd.concat([casedf, interpolatedf], ignore_index = True) #纵向拼接interpolateresult = pd.merge(data, result, how = "left") #匹配interpolateresult["diff"] = np.abs(interpolateresult["z"] - interpolateresult["lastmonthprice"])/interpolateresult["lastmonthprice"] #绝对偏差率for i in range(len(interpolateresult)): #利用插值更新本月基价if interpolateresult["label"][i] == 1 and interpolateresult["diff"][i]< 0.05: #如果插值的绝对误差率小于0.05填补空缺价格interpolateresult["thismonthprice"][i] = interpolateresult["z"][i]interpolateresult.to_excel("无案例楼盘interpolateresult.xlsx")return casenum, interpolatenum, interpolateresult[["newdiskid","plate", "newdiskname", "villa", "lon", "lat", "lastmonthprice", "thismonthprice"]] #返回更新本月基价的数据框if __name__ == "__main__":data = pd.read_excel("无案例楼盘逐步插值.xlsx") #读取数据scatterPlot(data)casenum, interpolatenum, data = gridInterpolate(data)
代码解读
主要由scatterPlot和gridInterpolate两个函数实现,
其中scatterPlot主要负责绘制散点图,图片如下
gridInterpolate主要负责插值功能
griddata基本调用格式如下
scipy.interpolate.griddata(points,values,xi,method ='linear',fill_value = nan,rescale = False )
- points
数据点坐标。可以是形状(n,D)的数组,也可以是ndim数组的元组。 - values
浮点或复数的ndarray,形状(n,)的数据值 - xi
浮点数的二维数组或一维数组的元组,形状(M,D)插值数据的点。
- method
可选插值方法。{‘linear’,‘nearest’,‘cubic’},之一,其中
linear 将输入点设置为n维单纯形,并在每个单形上线性插值,可以简单理解为以三角形为基础,就是按Delaunay方法先找出内插点四周的3个点,构成三角形,内插点在三角形内.然后线性内插,或三次方程内插.。
nearest 返回最接近插值点的数据点的值。
cubic 返回由三次样条确定的值。 返回由分段立方,连续可微(C1)和近似曲率最小化多项式表面确定的值。
- fill_value
float,可选,用于填充输入点凸包外部的请求点的值,如果未提供,则默认为nan。此选项对“最近”方法无效。 - rescale
bool,可选,在执行插值之前,重新缩放指向单位立方体,如果某些输入维度具有不可比较的单位并且相差很多个数量级。
逐步插值
插值是一个逐步扩散的过程,如果让第一次插值的结果再参与训练的话,第二次插值效果会好一些,以此类推,循环下去,就可以逐步插值,最后会达到一种收敛状态,所以需要用一个标志其达到收敛了,最简单的判断方式就是插值数据不再提升了就认为收敛了。
# -*- coding: utf-8 -*-
"""
Project_name:逐步插值
Description:多次插值实现无案例楼盘基价赋值
Created on Tue Oct 20 17:59:08 2020
@author: 帅帅de三叔
"""
import numpy as np #导入数值计算模块
import pandas as pd #导入数据分析模块
from scipy import stats
import matplotlib.pyplot as plt #导入绘图模块
from mpl_toolkits.mplot3d import Axes3D #导入3d绘图模块
from matplotlib import cm #颜色调整用
from scipy.interpolate import griddata #导入插值类def scatterPlot(data): #绘制绘图函数plt.figure(figsize = (9, 7)) #新建画布caseData = data.dropna(subset = ["thismonthprice"]) #去掉空值,保留本月有价格的数据留作插值依据lon = caseData["lon"] #经度lat = caseData["lat"] #纬度interpolateData = data[data["thismonthprice"].isnull()] #待插值的数据x = interpolateData["lon"] #待插值的楼盘经度y = interpolateData["lat"] #待插值的楼盘纬度plt.scatter(lon, lat, s = 5, color = "r", label = "no need interpolate") #不需要插值的散点图plt.scatter(x, y, s = 5, color = "b", alpha = 0.2, label = "need interpolate") #需要插值的散点图plt.legend()plt.show()def gridInterpolate(data): #定义插值函数,data需要有newdiskid, lon,lat, 上个月价格,本月价格(有空缺)casenum = [] #用来存放有价格的楼盘数interpolatenum = [] #用来存放待插值的样本数caseData = data.dropna(subset = ["thismonthprice"]) #去掉空值,保留本月有价格的数据留作插值依据casenum.append(len(caseData)) #追加有价格的样本数casediskid = caseData["newdiskid"] #楼盘idlon = caseData["lon"] #经度lat = caseData["lat"] #纬度z = caseData["thismonthprice"] #价格interpolateData = data[data["thismonthprice"].isnull()] #待插值的数据interpolatenum.append(len(interpolateData)) #追加待插值样本数interpolatediskid = interpolateData["newdiskid"] #待插值的楼盘idx = interpolateData["lon"] #待插值的楼盘经度y = interpolateData["lat"] #待插值的楼盘纬度xx, yy = np.meshgrid(x, y) #平面网格点待插值的位置zz = griddata((lon, lat), z, (xx, yy), method='nearest', fill_value = z.median(), rescale = True) #stats.mode(z)[0][0]在没有价格的网格点插值Z = pd.DataFrame(zz) #数据框化interpolateprice = [] #用来存放插值for i in range(Z.shape[0]):interpolateprice.append(Z[i][i]) #取对角线上的元素label1 = pd.Series(np.zeros(len(casediskid))) #有价格的标签casedf = pd.concat([casediskid, z], axis = 1).reset_index()[["newdiskid", "thismonthprice"]] #索引重置casedf = pd.concat([casedf, label1], axis = 1)casedf.columns = ["newdiskid", "z", "label"]label2 = pd.Series(np.ones(len(interpolatediskid))) #插值标签interpolatedf = pd.DataFrame({"interpolatediskid":interpolatediskid, "interpolateprice": interpolateprice}).reset_index()[["interpolatediskid", "interpolateprice"]]interpolatedf = pd.concat([interpolatedf, label2], axis = 1)interpolatedf.columns = ["newdiskid", "z", "label"] result = pd.concat([casedf, interpolatedf], ignore_index = True) #纵向拼接interpolateresult = pd.merge(data, result, how = "left") #匹配interpolateresult["diff"] = np.abs(interpolateresult["z"] - interpolateresult["lastmonthprice"])/interpolateresult["lastmonthprice"] #绝对偏差率for i in range(len(interpolateresult)): #利用插值更新本月基价if interpolateresult["label"][i] == 1 and interpolateresult["diff"][i]< 0.05: #如果插值的绝对误差率小于0.05填补空缺价格interpolateresult["thismonthprice"][i] = interpolateresult["z"][i]interpolateresult.to_excel("无案例楼盘interpolateresult.xlsx")return casenum, interpolatenum, interpolateresult[["newdiskid","plate", "newdiskname", "villa", "lon", "lat", "lastmonthprice", "thismonthprice"]] #返回更新本月基价的数据框if __name__ == "__main__":data = pd.read_excel("无案例楼盘逐步插值.xlsx") #读取数据scatterPlot(data)casenum, interpolatenum, data = gridInterpolate(data)print(casenum[-1], len(data.dropna(subset = ["thismonthprice"])))for i in range(10):if casenum[-1] < len(data.dropna(subset = ["thismonthprice"])): #未收敛则继续插值casenum, interpolatenum, data = gridInterpolate(data)scatterPlot(data)print(casenum[-1], len(data.dropna(subset = ["thismonthprice"])))
插值动态图如下,从动图可以看到需要插值的越来越少,即图中的蓝色的越来越少,反观红色的越来越多。
总结
插值还是很消耗资源的,所以比较慢,用3000个插值10000个,有种四两拨千斤的感觉,最后效果不会特别好,比如插值到收敛了一共插值出6000个,相当于填补了一部分数据标签,这为以后的机器学习模型提供了基础。
参考文献
1, 空间坐标和坐标所对应的属性(高程,温度等 )https://blog.csdn.net/csubai07/article/details/104344291
2, griddata用法 https://blog.csdn.net/gxf0789/article/details/82384372
3, griddata 原理 http://blog.sina.com.cn/s/blog_6264e23a01016k9f.html
4, matlab 的griddata https://www.mathworks.com/help/matlab/ref/griddata.html
5, Scipy 的插值 http://liao.cpython.org/scipytutorial11/
6,(数值分析)各种插值法的python实现 https://blog.csdn.net/qq_20011607/article/details/81412985
7,https://blog.csdn.net/weixin_43718675/article/details/103497930?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2allsobaiduend~default-1-103497930.nonecase&utm_term=python%E4%B8%ADgriddata%E5%87%BD%E6%95%B0&spm=1000.2123.3001.4430
griddata三维空间插值相关推荐
- 利用griddata进行插值
利用griddata进行插值 griddata函数讲解 第一步:导入相关库 第二步:给出插值到的经纬度信息(目标经纬度) 第三步:待插值数据 第四步:插值 汇总成函数 结果对比 插值前(10km) 插 ...
- python绘制气象等值线图_利用Python插值绘制等值线图
最近需要根据有限的站位点绘制插值等值线图,在网上中文搜索一通,只发现了这货Matplot Basemap 画湖北地图.插值.等值线,要么就是对这货的转载,这货不提供数据的形式,但是基本的代码思路还是不 ...
- python绘制等值线图_利用Python插值绘制等值线图
最近需要根据有限的站位点绘制插值等值线图,在网上中文搜索一通,只发现了这货Matplot Basemap 画湖北地图.插值.等值线,要么就是对这货的转载,这货不提供数据的形式,但是基本的代码思路还是不 ...
- 数学建模:插值与拟合—插值问题的python求解
目录 scipy.interpolate.interpnd 简介 一.一维插值 (interp1d) 二.二维网格节点插值 (interp2d) 三.二维散乱点插值 (griddata) scipy. ...
- python教程笔记(详细)
本文持续更新- 我的目录 说明 第一部分:python基础知识 0 python风格 0.1 命名风格 0.2 空行与空格 0.2.1 空格和tab 0.2.2 类中的空行问题 1 字符串,strin ...
- 【Python基础】科学计算库Scipy简易入门
0.导语 Scipy是一个用于数学.科学.工程领域的常用软件包,可以处理插值.积分.优化.图像处理.常微分方程数值解的求解.信号处理等问题.它用于有效计算Numpy矩阵,使Numpy和Scipy协同工 ...
- 【Matlab】离散点拟合曲面
Matlab中可以使用interp函数和griddata函数来实现插值,从而得到拟合曲面 从曲面上采样 离散点拟合曲面 1.离散点采样 为了获取离散点,首先需要从某个特定的曲面上采样. 其中利用数学的 ...
- matlab多项式操作
MATLAB 7.0 采用行向量来表示多项式,将多项式的系数按降幂次序存放在行向量中. 实例:设计一个函数 poly2str(),实现把一个行向量表示的多项式转换为常见的字符串表 示的多项式, %po ...
- 【机器学习基础】Scipy(科学计算库) 手把手手把手
0.导语 Scipy是一个用于数学.科学.工程领域的常用软件包,可以处理插值.积分.优化.图像处理.常微分方程数值解的求解.信号处理等问题.它用于有效计算Numpy矩阵,使Numpy和Scipy协同工 ...
- matlab空间曲面拟合,【Matlab】离散点拟合曲面
Matlab中可以使用interp函数和griddata函数来实现插值,从而得到拟合曲面 从曲面上采样 离散点拟合曲面 1.离散点采样 为了获取离散点,首先需要从某个特定的曲面上采样. 其中利用数学的 ...
最新文章
- 12门课100分直博清华!这份成绩单冲上热搜,但学霸小伙也曾考过25分
- Oracle错误大全(目前自己所碰到的)
- python asyncio_Python 的异步 IO:Asyncio 简介
- java中domain什么意思_java解析URL中domain、端口和协议的两种方法
- linux中循环创建文件,linux-尝试创建一个文件以调用另一个文件进行循环搜索
- [2019.3.21]洛谷P3640 [APIO2013]出题人
- 未来人工智能应用体现出的核心技术有哪些?
- hadoop component summary
- 拓端tecdat|R语言状态空间模型:卡尔曼滤波器KFAS建模时间序列
- 第9批候选!高工智能汽车金球奖入围年度产品/方案公示
- VS2015编译OpenDDS
- Elasticsearch:管理悬空(dangling)索引
- 李四光预测的地震带及合肥地震分析
- 全面了解电商网站建设要点,看这一篇就够了
- python-pandapower电力系统短路电流计算(算例3:探索一天的最佳电网运行方式))
- php 应用宝支付,手游渠道应用宝接入总结
- 算法实现之宝石与石头
- 邀请试用 实景三维模型在线浏览及网页分享平台
- 一起赚美金:经典Niche站变现模式分析(1)
- 计算机网络——应用层之文件传送协议(FTP)