前言
前段时间接了一个小程序活，对方似乎是需要卷积数据做实验，然后需要利用多个矢量点提取3X3核和5X5核点栅格数据做实验，委托我写一个python模块去提取这个所有矢量点周围的9个或25个栅格高程点，虽然我也不太理解这个实验，但是毕竟也不是没办法实现，所以索性接了下来，锻炼一下GDAL的使用。最后处理得到的结果见下图，效果还是可以的。
思路
这里借用了大佬的思路，核心思想是将数据进行仿射变换，转换到图上坐标以后，根据二维数组空间的加减获得卷积核，再利用构建好的值提取函数进行卷积核值提取，当然这么一说似乎很简单，但实际操作起来还是遇到不少问题的。
以上大佬的代码只够将值获取，但未输出栅格或矢量，这里加入自己的思路，在获取一个卷积核的9个坐标的情况下，利用OGR创建9个点的矢量点图层，再基于这个图层提取栅格值，最后使用矢量转栅格，即可得到以上部分内容的效果。
有了思路以后，开始写代码（google）
主函数文件
主函数文件主要放大佬写的仿射变换以及卷积核等函数，作为其他py文件导入使用。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
from osgeo import gdal
from osgeo import ogr,osr
import numpy as np
import re
#提取所有站点的经纬度
def get_vector_coordinates(filepath):
driver = ogr.GetDriverByName('ESRI Shapefile')
#打开矢量
ds = driver.Open(filepath, 0)
if ds is None:
print('Could not open ' +'stations.shp')
layer = ds.GetLayer()
#获取要素及要素地理位置
xValues = []
yValues = []
feature = layer.GetNextFeature()
while feature:
geometry = feature.GetGeometryRef()
x = geometry.GetX()
y = geometry.GetY()
xValues.append(x)
yValues.append(y)
feature = layer.GetNextFeature()
coordinates_list = np.array([x for x in zip(xValues,yValues)])
return coordinates_list
#生成提取模板
def Set_GridDataTemplate(grid_size = (5,5), scale = 0.1):
# 默认生成5*5的格网
grid = [[None] * grid_size[1] for i in range(grid_size[0])]
start_coordinate = -(grid_size[1] - 1) / 2 * scale
for i in range(grid_size[0]):
for j in range(grid_size[1]):
grid[i][j] = (start_coordinate + j * scale,-start_coordinate - i * scale)
template = np.array(grid,dtype=object)
return template
#得到站点周围的经纬度
def get_grid_coordinates(vector_coordinates,grid_template):
grid_coordinates_list = []
for i in vector_coordinates:
lon_lat_grid = [[i] * grid_template.shape[1] for j in range(grid_template.shape[0])]
grid_coordinates = lon_lat_grid + grid_template
grid_coordinates_list.append(grid_coordinates)
return grid_coordinates_list
#提取栅格值
def Extract_value(raster_path,grid_coordinates):
ds = gdal.Open(raster_path)
if ds is None:
print('Could not open image')
#获取行列、波段
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
#获取仿射变换信息
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]
all_feature = []
for i in grid_coordinates:
feature = []
count = 0
lon = 0
lat = 0
for j in i:
count += 1
for z in range(len(j)):
if (z + 1) == (len(grid_coordinates[0]) + 1)/2 and count == (len(grid_coordinates[0]) + 1)/2:
lon = j[z][0]
lat = j[z][1]
x = j[z][0]
y = j[z][1]
#获取点位所在栅格的位置
xOffset = int((x-xOrigin)/pixelWidth)
yOffset = int((y-yOrigin)/pixelHeight)
#提取每个波段上对应的像元值
for n in range(bands):
band = ds.GetRasterBand(n+1)
data = band.ReadAsArray(xOffset, yOffset,1,1)
if data != None:
value = int(data[0,0])
feature.append(value)
else:
feature.append(-999)
feature.append(lon)
feature.append(lat)
all_feature.append(feature)
return all_feature
if __name__ == "__main__":
shapefile = r'.\shapefile\stations.shp'
rasterfile = r'.\raster\dem.tif'
points_vector = get_vector_coordinates(shapefile)
grid55 = Set_GridDataTemplate(grid_size = (5,5), scale = 0.1)
coordinates_values = get_grid_coordinates(points_vector,grid55) #获取站点周围所有像元经纬度
cal_raster = Extract_value(rasterfile,coordinates_values)
#cal_raster2 = Extract_value(rasterfile,points_vector)
print(cal_raster)
#print(coordinates_values)
list_2 = [str(x) for item in coordinates_values for x in item]
out = open(r"raster.txt", "a+")
for i in cal_raster:
out.write(str(i))
out.close()
|
上述代码中，三个函数都是直接用的大佬写的内容，脚本文件运行部分是我自己加的，主要是将坐标点输出到raster.txt文本，如果想看详细的输出内容，可以把print()的注释去掉，将改python文件保存为func_all.py，结束第一步工作。当然这部分文件的脚本运行代码其实只是用来进行测试，后续其实主要是将该文件作为包导入。
OGR创建坐标txt文件
osgeo的OGR包，提供了创建ESRI shapefile的驱动器，直接手撸代码（google）就行，这个资源挺多。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
|
#导入相关包
from dataclasses import replace
from os import remove
from pkgutil import iter_modules
from matplotlib.pyplot import text
from osgeo import gdal
from osgeo import ogr,osr
import numpy as np
import re
import func_all
raster_path = r'\raster\dem.tif' #dem文件相对或绝对路径
filepath = r'\shapefile\stations.shp' #站点文件路径
driver = ogr.GetDriverByName('ESRI Shapefile')
#打开矢量
ds1 = driver.Open(filepath, 0)
layer = ds1.GetLayer()
xValues = []
yValues = []
feature = layer.GetNextFeature()
while feature:
geometry = feature.GetGeometryRef()
x = geometry.GetX()
y = geometry.GetY()
xValues.append(x)
yValues.append(y)
feature = layer.GetNextFeature()
coordinates_list = np.array([x for x in zip(xValues,yValues)])
grid55 = func_all.Set_GridDataTemplate(grid_size = (3,3), scale = 0.1) #这里grid55是做5*5实验的名字，实际运行的格网是3*3，但我也懒得改回来了
grid_coordinates_list = []
for i in coordinates_list:
lon_lat_grid = [[i] * grid55.shape[1] for j in range(grid55.shape[0])]
grid_coordinates = lon_lat_grid + grid55
grid_coordinates_list.append(grid_coordinates)
print(str(grid_coordinates_list).replace('[','').replace(']','') + "\n")
out = open(r"\raster3.txt", "a+") #根据需求修改路径输出坐标路径
for i in grid_coordinates_list:
text1 = str(i).replace('[','').replace(']','\n')
out.write(str(text1))
out.close()
# 调用文件的 readline()方法
file2 = open('\\text2.txt', 'w', encoding='utf-8') # 根据需求修改路径输出坐标文本路径
for line in file1.readlines():
line = line.lstrip(' ').replace(' ',',')
if line == "\n":
line = line.strip("\n")
file2.write(line)
num = line.split(',')
file2.close()
|
得到的txt文本结果如下图所示，该结果没有进行小数位保留，在后续操作中需要进行删位。
根据txt坐标文本创建矢量点图层
将该文件命名为create_point_shp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
|
from tokenize import Double
from numpy import double
from osgeo import ogr,osr,gdal
raster = r"\raster\dem.tif"
shape = r"\Point.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource(shape) ## shp文件名称
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) ## 空间参考：WGS84
layer = data_source.CreateLayer("Point", srs, ogr.wkbPoint) ## 图层名称要与shp名称一致
field_name = ogr.FieldDefn("Name", ogr.OFTString) ## 设置属性
field_name.SetWidth(20) ## 设置长度
layer.CreateField(field_name) ## 创建字段
field_Longitude = ogr.FieldDefn("Longitude", ogr.OFTReal) ## 设置属性
layer.CreateField(field_Longitude) ## 创建字段
field_Latitude = ogr.FieldDefn("Latitude", ogr.OFTReal) ## 设置属性
layer.CreateField(field_Latitude) ## 创建字段
field_height = ogr.FieldDefn("Height", ogr.OFTReal) ## 设置属性
layer.CreateField(field_height) ## 创建字段
#打开栅格
dr = gdal.Open(raster)
transform = dr.GetGeoTransform()
x_origin = transform[0]
y_origin = transform[3]
pixel_width = transform[1]
pixel_height = transform[5]
file1 = open('\\text2.txt', 'r', encoding='utf-8')
for line in file1.readlines():
currentline = line.split(',')
Longitude = float(currentline[0])
Latitude = float(currentline[1])
Longitude = round(Longitude,5) #保留5位小数
Latitude = round(Latitude,5)
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField("Name", "point") ## 设置字段值
feature.SetField("Longitude", str(Longitude)) ## 设置字段值
feature.SetField("Latitude", str(Latitude)) ## 设置字段值
x = Longitude
y = Latitude
x_offset = int((x - x_origin) / pixel_width)
y_offset = int((y - y_origin) / pixel_height)
dt = dr.ReadAsArray(x_offset, y_offset, 1, 1)
feature.SetField("Height", int(dt)) #添加高程列
wkt = "POINT(%f %f)" % (float(Longitude), float(Latitude)) ## 创建点
point = ogr.CreateGeometryFromWkt(wkt) ## 生成点
feature.SetGeometry(point) ## 设置点
layer.CreateFeature(feature) ## 添加点
feature.Destroy() ## 关闭属性
data_source.Destroy() ## 关闭数据
|
主要思路是，利用OGR的创建矢量数据的attributes，创建经纬列和高程值列，通过readline函数逐行读取坐标值，创建坐标点，再将根据该坐标点所对应的栅格数据高程值写入shp文件，即可实现根据坐标点的高程值提取至shp。效果见下图:
矢量转栅格
其实我认为，研究点周围的高程信息，直接用shp就够了，但客户要求要栅格，可能是想做别的实验吧，那就矢量转栅格。这部分再简单不过了，毕竟有手（google）就行，直接给代码，封装成了函数，直接用即可。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
|
from osgeo import gdal,ogr,gdalconst
def Rasterize(input_shp,input_tif,output_tif,field,filed_type,NoValue=-9999):
"""
input_shp:需要转为栅格的矢量文件（矢量文件路径）
input_tif:模板栅格，用于读取地理变换信息、栅格大小，将其应用于新的栅格上
output_tif:输出栅格文件（栅格文件路径）
field:字符串，栅格值的字段
filed_type:栅格值类型，一般选择gdal.GDT_Int16,gdal.GDT_Int32,gdal.GDT_Float32,gdal.GDT_Float64等几种类型
NoValue:整型或浮点型，矢量空白区转换后的值
"""
data = gdal.Open(input_tif, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
proj=data.GetProjection()
open_shp = ogr.Open(input_shp)
shp_ly = open_shp.GetLayer()
x_min, x_max, y_min, y_max = shp_ly.GetExtent()
pixel_size = geo_transform[1]
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(output_tif, x_res, y_res, 1, filed_type)
target_ds.SetGeoTransform((x_min, pixel_size, 0.0, y_max, 0.0, -pixel_size))
target_ds.SetProjection(proj)
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoValue)
band.FlushCache()
if field is None:
gdal.RasterizeLayer(target_ds, [1], shp_ly, None)
else:
OPTIONS = ['ATTRIBUTE=' + field]
gdal.RasterizeLayer(target_ds, [1], shp_ly, options=OPTIONS)
target_ds = None
input_tif1 = r"\raster\dem.tif"
input_shape = r"\Point.shp"
output_tif1 = r"\test.tif"
Rasterize(input_shp=input_shape,input_tif=input_tif1,output_tif=output_tif1,field= 'Height',filed_type=gdal.GDT_Int16,NoValue=-9999)
|
最终得到一个全国分布的卷积核3X3图像，效果见下图，细节效果看第一节内容的图。
后记
如当然这个项目还有很多不完善的地方，比如其实用到的函数都放在开头说的func_all比较好，但小项目嘛，把代码给到客户，客户说还需要自己改造，也就算了。
实施的思路其实还是不够简洁，目前想到如果在提取坐标以后，不输出到txt，而直接用正则表达式re衔接到第二部分，是不是会提高效率一些？
如果在提取栅格到矢量点的过程，直接利用坐标值和高程值创建新栅格，是不是更快一些？这样就少了过程的shp数据，当然目前我还不知道如何直接利用这些信息构建栅格，思路没打开还是难😅😅
另外说一下，这里的实验所有使用的数据，全都是来自大佬的网盘，密码是dkbw。
参考博文
- 从WKT创建几何
- 创建新的shp数据并添加数据
- 栅格化矢量图层
- Python GDAL 地学分析（八）提取矢量点周围栅格值
