GDAL实现核密度分析

from math import log,sqrt,pi
import geopy
from osgeo import gdal
from osgeo import ogr
import numpy as np
from osgeo import osr
from osgeo.gdalconst import GDT_Float32
from geopy import distance
import pyproj


class KernelDensity(object):
    def __init__(self,inputPath:str,outPath:str,cellSize=None,population=None,searchRadius=None) -> None:
        super().__init__()
        self.inputPath=inputPath
        self.outPath=outPath
        self.population=population
        self.cellSize=cellSize
        self.radius=searchRadius
        self.__read()
        self.__radius()
        self.__density()
        self.__outTif()
        print('核密度计算完毕！')
    def __read(self):
        ogr.RegisterAll()
        ds=ogr.Open(self.inputPath)
        self.ds=ds
        lyr=ogr.DataSource.GetLayer(self.ds,0)
        xList=list()
        yList=list()
        pList=list()
        if self.population:
            index=ogr.Layer.FindFieldIndex(lyr,self.population,True)
        else:
            index=-1
        while True:
            feature=ogr.Layer.GetNextFeature(lyr)
            if not(feature):
                break
            geometry=ogr.Feature.geometry(feature)
            x=ogr.Geometry.GetX(geometry)
            y=ogr.Geometry.GetY(geometry)
            xList.append(x)
            yList.append(y)
            if index==-1:
                pList.append(1)
            else:
                pList.append(ogr.Feature.GetField(feature,index))
        xList=np.array(xList)
        yList=np.array(yList)
        pList=np.array(pList)
        pSum=np.sum(pList)
        wList=pList/pSum
        self.xList=xList
        self.yList=yList
        self.wList=wList
        self.pList=pList
    def __radius(self):
        # 计算默认搜索半径
        if not(self.radius):
            xC=np.dot(self.xList,self.wList)
            yC=np.dot(self.yList,self.wList)
            lyr=ogr.DataSource.GetLayer(self.ds,0)
            sr=ogr.Layer.GetSpatialRef(lyr)
            xyList=zip(self.xList,self.yList)
            dList=[KernelDensity.dist([xC,yC],xy,sr) for xy in xyList]
            dList=np.array(dList)
            SD=sqrt(np.dot(dList**2,self.wList))   
            # 计算要素点到中心点的距离同样需要加权
            dist=dList*self.pList
            distS=np.sort(dist)
            count=distS.shape[0]
            if count%2:
                Dm=distS[count//2]
            else:
                Dm=(distS[count//2]+distS[count//2-1])/2
            radius=min(SD,sqrt(1/log(2))*Dm)
            if self.population:
                radius=0.9*radius*pow(np.sum(self.pList),-0.2)
            else:
                radius=0.9*radius*pow(count,-0.2)
            self.radius=radius
            self.sr=sr
        print("默认搜索半径为:{}".format(radius))
        return
    @staticmethod
    def dist(p1,p2,sr:osr.SpatialReference):
        x1,y1=p1
        x2,y2=p2
        if sr.IsProjected():
            dx=x1-x2
            dy=y1-y2
            return sqrt(dx**2+dy**2)
        else:
            # 将地理坐标统一到wgs84下便于计算
            srProj=osr.SpatialReference.ExportToProj4(sr)
            sourceProj=pyproj.Proj(srProj)
            dstProj=pyproj.Proj('epsg:4326')
            x1,y1=pyproj.transform(sourceProj,dstProj,x1,y1)
            x2,y2=pyproj.transform(sourceProj,dstProj,x2,y2)
            dist=distance.geodesic((x1,y1),(x2,y2)).meters
            return dist
    def __density(self):
        lyr=ogr.DataSource.GetLayer(self.ds,0)
        lyrExtent=ogr.Layer.GetExtent(lyr)
        width=lyrExtent[1]-lyrExtent[0]
        height=lyrExtent[3]-lyrExtent[2]
        if not(self.cellSize):
            self.cellSize=min(width,height)/250
        xSize=int(width/self.cellSize)
        ySize=int(height/self.cellSize)
        densityArray=np.zeros((ySize,xSize))
        yMax=lyrExtent[3]
        xMin=lyrExtent[0]
        self.topleft=[xMin,yMax]
        # num 输入数据集中的点数
        num=ogr.Layer.GetFeatureCount(lyr)
        for i in range(ySize):
            y=yMax-i*self.cellSize-0.5*self.cellSize
            for j in range(xSize):
                x=xMin+j*self.cellSize+0.5*self.cellSize
                density=0
                for k in range(num):
                    # 数据集坐标为投影坐标
                    xP=self.xList[k]
                    yP=self.yList[k]
                    dx=xP-x
                    dy=yP-y
                    dist=KernelDensity.dist([x,y],[xP,yP],self.sr)
                    if dist<self.radius:
                        density+=3/pi*self.pList[k]*pow((1-pow(dist/self.radius,2)),2)
                density*=pow(self.radius,-2)
                densityArray[i,j]=density*pow(10,6)
        self.densityArray=densityArray
        return
    def __outTif(self):
        # 将核密度计算的结果输出
        gdal.AllRegister()
        tifDriver=gdal.GetDriverByName('GTiff')
        row,col=self.densityArray.shape
        outds=gdal.Driver.Create(tifDriver,self.outPath,col,row,1,GDT_Float32,['COMPRESS=LZW'])
        outband=gdal.Dataset.GetRasterBand(outds,1)
        gdal.Band.WriteArray(outband,self.densityArray,0,0)
        # 设置波段的坐标参考
        lyr=ogr.DataSource.GetLayer(self.ds,0)
        crs=ogr.Layer.GetSpatialRef(lyr)
        gdal.Dataset.SetSpatialRef(outds,crs)
        geotransform=[self.topleft[0],self.cellSize,0,self.topleft[1],0,-self.cellSize]
        gdal.Dataset.SetGeoTransform(outds,geotransform)
        return
kd=KernelDensity('./Supermarket.shp','out6.tif',10,'money',None)
# kd=KernelDensity('./Supermarket.shp','out5.tif',10,None,None)`!

在这里插入图片描述