gdal 实现 热点分析 (Getis-Ord Gi*)

算法实现参考ArcGIS 对 Getis-Ord Gi* 的介绍
http://desktop.arcgis.com/zh-cn/arcmap/latest/tools/spatial-statistics-toolbox/h-how-hot-spot-analysis-getis-ord-gi-spatial-stati.htm

#include "TDGetisOrdGi.h"
#include "ogrsf_frmts.h"
#include <algorithm>
#include "TDAPICommon.h"
bool TDGetisOrdGi::process(LayerPathInfoV1 layerpathinfo, LayerPathInfoV1 armlayerpathinfo, std::string tempFields)
{
	m_selectFields = tempFields;
	auto fun = std::bind(&TDGetisOrdGi::cal, this);
	return process_byfunction(layerpathinfo, armlayerpathinfo, fun);
}
void TDGetisOrdGi::setLayerFeatureDfn(){
	m_ogrfeaturedefn->AddFieldDefn(new OGRFieldDefn("Gi", OFTReal));
}
void TDGetisOrdGi::cal()
{
	m_orgLayer->ResetReading();
	OGRFeature *poFeature;
	while ((poFeature = m_orgLayer->GetNextFeature()) != NULL)
	{
		OGRIntermediatedata tempogrintermediatedata;
		std::get<0>(tempogrintermediatedata) = poFeature->GetFID();
		std::get<1>(tempogrintermediatedata) = poFeature->GetFieldAsDouble(m_selectFields.c_str());
		auto pgeom = poFeature->GetGeometryRef();
		if (NULL != pgeom)
		{
			pgeom->Centroid(&std::get<2>(tempogrintermediatedata));
		}
		m_OGRIntermediatedata.push_back(tempogrintermediatedata);
	}
	std::for_each(m_OGRIntermediatedata.begin(), m_OGRIntermediatedata.end(), [&](OGRIntermediatedata & obj) {
		m_mapG[std::get<0>(obj)] = calGetisOrdGi(obj);
	});
	m_orgLayer->ResetReading();
	while ((poFeature = m_orgLayer->GetNextFeature()) != NULL)
	{
		auto NewFeature = poFeature->Clone();
		NewFeature->SetField("Gi", m_mapG[poFeature->GetFID()]);
		m_cachevector.push_back(NewFeature);
	}
	wiriteInfo();
}
double TDGetisOrdGi::calGetisOrdGi(OGRIntermediatedata obj)
{
	//求wij * xj 和
	double Gi = 0.0;
	double sumwijxj = 0.0;
	double sumwij = 0.0;
	double sumxj = 0.0;
	double sumwij2 = 0.0;
	double sumxj2 = 0.0;
	int n = m_OGRIntermediatedata.size();
	std::for_each(m_OGRIntermediatedata.begin(), m_OGRIntermediatedata.end(), [&](OGRIntermediatedata & tempobj) {
		if (std::get<0>(tempobj) != std::get<0>(obj))
		{
			double wij = TDAPICommon::computeDistanceBearing(std::get<2>(obj), std::get<2>(tempobj), pmoSourceSRS);
			double xj = std::get<1>(tempobj);
			sumwijxj += wij*xj;
			sumwij += wij;
			sumxj += xj;
			sumwij2 += std::pow(wij,2);
			sumxj2 += std::pow(xj, 2);
		}
	});
	double X = sumxj / n;
	double S = std::sqrt((sumxj2/n) - std::pow(X,2));
	Gi = (sumwijxj - (X*sumwij)) / (S*std::sqrt((n*sumwij2-std::pow(sumwij,2))/n-1));
	return Gi;
}

基类可以参考前几篇博客，算法实现是通过读取某列值，获取模型中心点将中心点和FID存下来，通过calGetisOrdGi计算出来每个要素的Gi*值，通过标准正太分布来圈定高低聚类
结果:

计算结果上来看还是挺靠谱的