【Python&GIS】通过经纬度创建矢量点文件
最近在做项目时，需要判断某个点是否在感兴趣区内。所以需要使用Python先根据经纬度的点创建矢量文件，再通过点文件和面文件的位置关系判断点是否在面内。
这里我们使用osgeo中的ogr和osr库，ogr库是一个处理地理空间矢量数据的开源库。它可以读取多种数据格式，进行地理处理、属性表操作、数据分析等操作。目前ogr和osr库已集成到GDAL库中，可以对栅格数据、矢量数据进行处理分析，被3S的研究人员广泛应用。感兴趣的可以自己去了解一下，不懂得可以一起交流！
1.安装所需库包，因为GDAL、OGR、OSR以及合并成osgeo中，所以需要从osgeo中导入。
from osgeo import ogr, osr2.初始化资源空间，这里需要写入保存的目标目录及文件名
driver = ogr.GetDriverByName("ESRI Shapefile")
# 创建数据驱动
source_data = driver.CreateDataSource("G:/Point1.shp")
# 创建数据资源3.创建投影空间，设置投影信息。这里我的点是WGS84的地理坐标系，所以就设置输出的文件为WGS84地理坐标系。（二选一即可）
1）4386是WGS84的EPSG编码，库里面内置了很多坐标系的EPSG编码。可以在库包中找到一些EPSG编码。存放在“C:\Program Files\Python36\Lib\site-packages\osgeo\data\gdal”中的ozi_datum.csv文件中，其他的可以自己百度搜索，也可以去EPSG官网查询。
spatial_proj = osr.SpatialReference()
# 创建SpatialReference对象，再向SpatialReference导入投影信息
spatial_proj.ImportFromEPSG(4326)2）当然我们也可以从已有的坐标系中导入，例如我们可以导入一个包含坐标系的矢量文件，读取它的坐标系信息，使用这个坐标系作为输出的坐标系（这一步可以导入一些地方坐标系，当然你需要有包含这个坐标系的矢量文件）。
file_path = "G:/local.shp"
ds = ogr.Open(file_path) # 打开数据集dataset
layer_ds = ds.GetLayer()
proj_ds = layer_ds.GetSpatialRef()
# 读取已有坐标系4.创建图层，创建特征信息（属性字段），编辑字段名设置字段长度等。
layer = source_data.CreateLayer("point", spatial_proj, ogr.wkbPoint)
# 创建图层，保持名称与文件名一致
field_longitude = ogr.FieldDefn("Longitude", ogr.OFTReal)
field_latitude = ogr.FieldDefn("Latitude", ogr.OFTReal)
# 创建字段，文本属性
field_longitude.SetWidth(10)
field_latitude.SetWidth(10)
# 设置字段长度
layer.CreateField(field_longitude)
layer.CreateField(field_latitude)
# 创建字段5.写入字段信息，创建点的几何位置。记得要关闭属性表和数据资源，不然数据在内存中不会保存。如果有多个点，可以创建列表保存经纬度，使用for循环遍历写入属性表。
feature_point = ogr.Feature(layer.GetLayerDefn())
# 创建feature
feature_point.SetField("Longitude", 126.123)
feature_point.SetField("Latitude", 31.123)
# 输入字段值
point_geo = ogr.Geometry(ogr.wkbPoint)
# 创建几何点
point_geo.AddPoint(126.123, 31.123)
# 添加几何点
feature_point.SetGeometry(point_geo)
# 设置点的字段值
layer.CreateFeature(feature_point)
feature_point.Destroy()
# 关闭属性
source_data.Destroy()6. 完整代码
根据自己的情况修改坐标系统以及文件保存的目录。第三步的二选一，我默认使用的是内置的EPSG编码，导入已有的投影信息已经注释掉了。如果需要使用已有的投影信息，将内置的EPSG两行代码删掉，再删掉注释的"就行了。
我这里是有多个点需要创建，所以已经使用了for循环去遍历经纬度列表，根据需求自行修改。
# -*- coding: utf-8 -*-
"""
@Time ： 2023/5/19 9:05
@Auth ： RS迷途小书童
@File ：Create Multipoint.py
@IDE ：PyCharm
"""
from osgeo import ogr, osr
from pyproj import Proj
def Create_multipoint(list_longitude, list_latitude, path_result):
"""
:param list_longitude: 输入经度列表
:param list_latitude: 输入纬度列表
:param path_result: 输入保存的shp路径
:return: 返回shp路径
"""
driver = ogr.GetDriverByName("ESRI Shapefile")
# 创建数据驱动
source_data = driver.CreateDataSource(path_result)
# 创建数据资源
spatial_proj = osr.SpatialReference()
# 创建SpatialReference对象，再向SpatialReference导入投影信息
spatial_proj.ImportFromEPSG(4326)
# EPSG编码：【4326：WGS84 、 32651：UTM/WGS84 51N】
"""file_path = "G:/local.shp"
ds = ogr.Open(file_path) # 打开数据集dataset
layer_ds = ds.GetLayer()
proj_ds = layer_ds.GetSpatialRef()
# 读取已有坐标系"""
layer = source_data.CreateLayer("point", spatial_proj, ogr.wkbPoint)
# 创建图层，保持名称与文件名一致
field_longitude = ogr.FieldDefn("Longitude", ogr.OFTReal)
field_latitude = ogr.FieldDefn("Latitude", ogr.OFTReal)
# 创建字段，文本属性
field_longitude.SetWidth(10)
field_latitude.SetWidth(10)
# 设置字段长度
layer.CreateField(field_longitude)
layer.CreateField(field_latitude)
# 创建字段
for i in range(len(list_longitude)):
feature_point = ogr.Feature(layer.GetLayerDefn())
# 创建feature
feature_point.SetField("Longitude", str(list_longitude[i]))
feature_point.SetField("Latitude", str(list_latitude[i]))
# 输入字段值
point_geo = ogr.Geometry(ogr.wkbPoint)
# 创建几何点
point_geo.AddPoint(float(list_longitude[i]), float(list_latitude[i]))
# 添加几何点
feature_point.SetGeometry(point_geo)
# 设置点的字段值
layer.CreateFeature(feature_point)
feature_point.Destroy()
# 关闭属性
source_data.Destroy()
# 关闭数据
return path_result
if __name__ == "__main__":
list_longitude = [100.1, 100.2, 100.3, 100.4]
list_latitude = [30.1, 30.2, 30.3, 30.4]
path_result = "G:/Point1.shp"
Create_multipoint(list_longitude, list_latitude, path_result)
本文章主要是分享个人在学习Python过程中写过的一些代码。有些部分借鉴了前人以及官网的教程，如有侵权请联系作者删除，大家有问题可以随时留言交流，博主会及时回复。
