python实现对csv文件对地理数据进行核密度估计

网上关于核密度估计的算法大多是一维数据，而基于地理数据的核密度估计相对较少，这里查阅了网上一些资料，自己改了一些地方，对带经纬度的csv数据进行了核密度估计。

首先先加载需要的包，开源的python包当中肯定必不可少的就是gdal库，这里除了gdal库还包括geopandas，目的主要是对csv文件转成矢量数据相对简洁一点点；pyproj库主要是用于投影转换的，而math库是用于数学计算时对一些数学符号的调用。

import pandas as pd
import geopandas
import pyproj
from osgeo import gdal, gdalconst
from osgeo import ogr
import numpy as np
import math

 那么开始了，第一步，导入csv文件，并进行投影转换，最后得到shp数据。

df= geopandas.read_file('E:/zjg_data/battery2015.csv',encoding='utf-8')#加载csv文件
df[['Longitude','Latitude']] = df[['Longitude','Latitude']].apply(pd.to_numeric)
gdf = geopandas.GeoDataFrame(df, geometry=geopandas.points_from_xy(df.Longitude, df.Latitude)) #将csv文件中的两个字段设为经纬度，转换成点
gdf.insert(4, 'val', 1, allow_duplicates=False)#插入一个字段val，值为1
gdf.crs = pyproj.CRS.from_user_input('EPSG:4326') #给输出的shp定义投影
gdf.to_crs(epsg=32617, inplace=True)#投影转换
gdf.to_file('E:/zjg_data/try/py_data/111/outpoint.shp', driver='ESRI Shapefile',encoding='utf-8')#输出shp数据

 紧接着，对上一步计算得到的shp数据进行栅格化，注意一下这里的rasterFile指定是csv数据所在的区域的栅格数据。目的是为了为shp转栅格数据提供栅格数据的一些参数。最后输出得到的是栅格点

shpFile='E:/zjg_data/try/py_data/111/outpoint.shp' #点矢量数据
rasterFile = 'E:/zjg_data/try/inareapro2_1.tif'  # 原影像

dataset = gdal.Open(rasterFile, gdalconst.GA_ReadOnly)#打开栅格数据

geo_transform = dataset.GetGeoTransform()#利用GetGeoTransfrom函数可以获得图像中任意一点的地理坐标

cols = dataset.RasterXSize  # 列数
rows = dataset.RasterYSize  # 行数
#影像左上角横坐标：geoTransform[0]，影像左上角纵坐标：geoTransform[3]遥感图像的水平空间分辨率为geoTransform[1]，遥感图像的垂直空间分辨率为geoTransform[5]，通常geoTransform[5] 与 geoTransform[1]相等
#如果遥感影像方向没有发生旋转，即上北、下南，则geoTransform[2] 与 row *geoTransform[4] 为零。
x_min = geo_transform[0]
y_min = geo_transform[3]
pixel_width = geo_transform[1]

shp = ogr.Open(shpFile, 0)
m_layer = shp.GetLayerByIndex(0)#获得第一个图层
target_ds = gdal.GetDriverByName('GTiff').Create("E:/zjg_data/try/py_data/111/raster_point.tif", xsize=cols, ysize=rows, bands=1, eType=gdal.GDT_Byte)#新创一个栅格文件
#对新创的栅格数据设置和原来栅格数据相同的地理坐标和投影
target_ds.SetGeoTransform(geo_transform)
target_ds.SetProjection(dataset.GetProjection())

band = target_ds.GetRasterBand(1)#获得第一个波段
band.SetNoDataValue(0)#设置无效值
band.FlushCache()#清理内存

gdal.RasterizeLayer(target_ds, [1], m_layer, options=["ATTRIBUTE=val"])  # 跟shp字段给栅格像元赋值
# gdal.RasterizeLayer(target_ds, [1], m_layer) # 多边形内像元值的全是255
del dataset
del target_ds

下一步，统计需要进行核密度估计的点的数目和位置以及设置核密度估计的搜索半径，并设置一个空的矩阵，等待数据的计算。

dataset2 = gdal.Open("E:/zjg_data/try/py_data/111/raster_point.tif", gdalconst.GA_ReadOnly)#打开栅格点数据
dataset2_geo_trans = dataset2.GetGeoTransform()
x_min2 = dataset2_geo_trans[0]
y_max2 = dataset2_geo_trans[3]
x_res2 = dataset2_geo_trans[1]
col = dataset2.RasterXSize
row = dataset2.RasterYSize

min_col_row = col if col < row else row#确定图像长和宽谁最大
radius=3*x_res2#设置搜索半径
#radius = (min_col_row / 30) * x_res2
radius_pixel_width = int(radius / abs(x_res2))#设置搜索半径的网格数

raster_data = dataset2.GetRasterBand(1).ReadAsArray()#以数组的形式读取栅格点数据
points_xy = []
values = []
point_shp_ds = ogr.GetDriverByName("ESRI Shapefile").Open(shpFile, 0)
if shpFile is None:
    print("文件{0}打开失败！".format(shpFile))

points_layer = point_shp_ds.GetLayer(0)
feature_count = points_layer.GetFeatureCount()#读取矢量点数据中的元素（点）数量
for i in range(feature_count):
    feature = points_layer.GetFeature(i)
    feature_defn = points_layer.GetLayerDefn()#获取数据表信息
    calc_field_index = feature_defn.GetFieldIndex("val")#val作为属性值
    field_value = feature.GetField(calc_field_index)
    values.append(field_value)

    point_xy = feature.GetGeometryRef()#提取feature的几何形状
    points_xy.append([point_xy.GetX(), point_xy.GetY()])#设置x，y

out_raster2="E:/zjg_data/try/py_data/111/outraster3.tif"
out_ds = dataset2.GetDriver().Create(out_raster2, col, row, 1,
                                        dataset2.GetRasterBand(1).DataType)  # 创建一个构建重采样影像的句柄
out_ds.SetProjection(dataset2.GetProjection())  # 设置投影信息
out_ds.SetGeoTransform(dataset2_geo_trans)  # 设置地理变换信息
result_data = np.full((row, col), -999.000, np.float32)  # 设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值

最后，进行核密度估计，这里采用的核函数是高斯函数

pi=math.pi#设置Π
print(len(points_xy))
for i in range(len(points_xy)):

    x_pos = int((points_xy[i][0] - x_min2) / abs(x_res2))
    y_pos = int((y_max2 - points_xy[i][1]) / abs(x_res2))
    # 先搜索以点为中心正方形区域，再在里面搜索圆形
    y_len_min = (y_pos - radius_pixel_width - 1) if (y_pos - radius_pixel_width - 1) > 0 else 0
    y_len_max = (y_pos + radius_pixel_width + 1) if (y_pos + radius_pixel_width + 1) < row else row

    x_len_min = (x_pos - radius_pixel_width - 1) if (x_pos - radius_pixel_width - 1) > 0 else 0
    x_len_max = (x_pos + radius_pixel_width + 1) if (x_pos + radius_pixel_width + 1) < col else col
    for y in range(y_len_min, y_len_max):
        for x in range(x_len_min, x_len_max):
            distance = ((x - x_pos) ** 2 + (y - y_pos) ** 2) ** 0.5
 				# 判断在半径内
            if (distance < radius_pixel_width and distance > 0):
                value = raster_data[y_pos][x_pos]
                k=math.exp(-(distance**2)/(2*radius_pixel_width))
                D_value=k/((radius_pixel_width**2)*pi*math.sqrt(2*pi))

                #D_value=(distance/radius_pixel_width) / ((radius_pixel_width**2)*pi)
                #if (raster_data[y][x] != 0):
                if (result_data[y][x] != -999.0):
                    result_data[y][x] += D_value
                else:
                    result_data[y][x] = D_value
            elif(distance==0):
                result_data[y][x]=raster_data[y][x]

result_data[result_data<0]=0
result_data=100*result_data
datamax=result_data.max()
datamin=result_data.min()
print(datamax,datamin)
dataResult=255*(result_data/datamax)
result_data[result_data==0]=-999

out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(result_data)
out_band.SetNoDataValue(-999)
out_band.FlushCache()
out_band.ComputeStatistics(False)  # 计算统计信息
out_ds.BuildOverviews('average', [1, 2, 4, 8, 16, 32])  # 构建金字塔

del out_ds  # 删除句柄
print("运行结束")

文章借鉴的文章：

使用python+gdal实现arcgis的核密度分析，针对点密度_gis_rc的博客-CSDN博客