Python|矿产卫片Excel经纬度坐标数据转换为shp点数据——OGR库实现

1.实验需求

基于Excel表格里面的经纬度坐标数据，自动生成点shp矢量文件，并添加属性信息。

2.编程思路详解

①使用Pandas库读取原始矿产图斑列表表格；

xlsx_path = u'C:\\Users\\YaoJun\\Desktop\\矿产图斑列表.xlsx'
#sheet_name默认为0，即读取第一个sheet的数据
df = pd.read_excel(xlsx_path, sheet_name=0, index_col=0, skiprows=2)

②将中心点坐标这一列拆分为X坐标和Y坐标两列，分别去除X:/Y:多余字符；

#将中心点坐标列拆分为两列
df[[u'X坐标', u'Y坐标']] = df[u'中心点坐标'].str.split(expand=True)
#去掉X:与Y：
df[u'X坐标'] = df[u'X坐标'].map(lambda x: x.replace('X:', ''))
df[u'Y坐标'] = df[u'Y坐标'].map(lambda x: x.replace('Y:', ''))

③分别将X坐标和Y坐标两列度分秒格式转换为十进制格式；

#度分秒转换为经纬度坐标
df[u'X坐标'] = df[u'X坐标'].apply(longitude_, )
df[u'Y坐标'] = df[u'Y坐标'].apply(longitude_, )

其中，度分秒转十进制坐标的公式如下：

十进制度数 = 度数 + 分数/60 + 秒数/3600。其中，度数为整数部分，分数为小数部分的整数部分，秒数为小数部分的小数部分。

例如：经度为120°30'15"，则十进制度数为：120 + 30/60 + 15/3600 = 120.50416666666667

度分秒转为十进制函数如下：

#度分秒转经纬度
def longitude_(longitude):
    longitude_split = re.split(u"°|\′|\″|\'|\"", longitude)[:3]
    if len(longitude_split) == 3:
        x = [float(j) for j in longitude_split]
        data = (x[0] + x[1] / 60 + x[2] / 3600)
        return str('%.6f'% float(data))

④使用OGR库创建与Excel文件同名的点shp文件，并将表格所有列名称添加为shp字段；

 #创建输出shp文件
 shpfname = xlsx_path.split('.')[0] + ".shp"
 CreatePolygonShp(shpfname)

创建shp函数：

#创建shp文件
def CreatePolygonShp(shpfname, epsg = 4490):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(shpfname):
        driver.DeleteDataSource(shpfname)
    ds = driver.CreateDataSource(shpfname)
    SpatialReference = osr.SpatialReference()
    SpatialReference.ImportFromEPSG(epsg)
    #ds.CreateLayer(shpfname, SpatialReference, geom_type=ogr.wkbPoint, options=['ENCODING=GBK'])
    ds.CreateLayer(shpfname, SpatialReference, geom_type=ogr.wkbPoint)
    del ds
    print("1.shp文件创建成功！")

创建字段函数：

#给指定shp文件添加字段函数
def AddField(shpfname, fieldname, fieldtype, fieldlenth = 50):
    ds = ogr.Open(shpfname, 1)
    layer = ds.GetLayer(0)
    # 添加字段
    fieldDefn = ogr.FieldDefn(fieldname, fieldtype)
    if fieldtype == ogr.OFTString:
        fieldDefn.SetWidth(fieldlenth)
    layer.CreateField(fieldDefn)
    del ds

将表格所有列名称添加为shp文件字段：

#添加所有列名为字段
    for col_name in df.columns[:-2]:
        if len(col_name) > 5:
            continue
        if u"面积" in col_name:
            AddField(shpfname, col_name, ogr.OFTReal)
        else:
            AddField(shpfname, col_name, ogr.OFTString, 30)

⑤遍历表格所有行，依次创建点图形并添加属性信息。

#打开空白shp
driver = ogr.GetDriverByName('ESRI Shapefile')
fc1 = driver.Open(shpfname, 1)
layer = fc1.GetLayer()
layer_defn = layer.GetLayerDefn()
field_count = layer_defn.GetFieldCount()

#添加图形与属性信息
for i in range(0, len(df.index)):
    feature = ogr.Feature(layer.GetLayerDefn())
    #依次给字段附属性值
    for j in range(field_count):
        field_defn = layer_defn.GetFieldDefn(j)
        field_name = field_defn.GetName()
        if field_name in df.columns:
            feature.SetField(field_name, str(df.iloc[i][field_name]))
    #创建点
    wkt = "POINT({} {})".format(df.iloc[i][u'X坐标'], df.iloc[i][u'Y坐标'])  # "point({} {})"两个参数之间没有逗号，记住
    point = ogr.CreateGeometryFromWkt(wkt)
    #print(point)
    feature.SetGeometry(point)
    layer.CreateFeature(feature)
    del feature

3.完整代码

# -*- coding: utf-8 -*-
# @Time    : 2023/4/16 15:49
# @Author  : 药菌

import re
import pandas as pd
from osgeo import ogr,osr,gdal
import os


gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")  # 支持文件名称及路径内的中文
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")  # 支持属性字段中的中文

#度分秒转经纬度
def longitude_(longitude):
    longitude_split = re.split(u"°|\′|\″|\'|\"", longitude)[:3]
    if len(longitude_split) == 3:
        x = [float(j) for j in longitude_split]
        data = (x[0] + x[1] / 60 + x[2] / 3600)
        return str('%.6f'% float(data))

#经纬度转度分秒
def Degrees(longitude):
    data1 = int(float(longitude))
    data2 = float(float('0.'+longitude.split('.')[1])*60)
    data3 = int(data2)
    data4 = float('0.'+str(data2).split('.')[1])*60
    data5 = '%.2f'% float(data4)
    data6 = int(float(data5))
    print('度',data1)
    print('分',data3)
    print('秒',data5)
    return [data1,data3,data5]

#创建shp文件
def CreatePolygonShp(shpfname, epsg = 4490):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(shpfname):
        driver.DeleteDataSource(shpfname)
    ds = driver.CreateDataSource(shpfname)
    SpatialReference = osr.SpatialReference()
    SpatialReference.ImportFromEPSG(epsg)
    #ds.CreateLayer(shpfname, SpatialReference, geom_type=ogr.wkbPoint, options=['ENCODING=GBK'])
    ds.CreateLayer(shpfname, SpatialReference, geom_type=ogr.wkbPoint)
    del ds
    print("1.shp文件创建成功！")

#给指定shp文件添加字段函数
def AddField(shpfname, fieldname, fieldtype, fieldlenth = 50):
    ds = ogr.Open(shpfname, 1)
    layer = ds.GetLayer(0)
    # 添加字段
    fieldDefn = ogr.FieldDefn(fieldname, fieldtype)
    if fieldtype == ogr.OFTString:
        fieldDefn.SetWidth(fieldlenth)
    layer.CreateField(fieldDefn)
    del ds

if __name__ == '__main__':
    xlsx_path = u'C:\\Users\\YaoJun\\Desktop\\矿产图斑列表.xlsx'
    #sheet_name默认为0，即读取第一个sheet的数据
    df = pd.read_excel(xlsx_path, sheet_name=0, index_col=0, skiprows=2)

    #将中心点坐标列拆分为两列
    df[[u'X坐标', u'Y坐标']] = df[u'中心点坐标'].str.split(expand=True)

    #去掉X:与Y：
    df[u'X坐标'] = df[u'X坐标'].map(lambda x: x.replace('X:', ''))
    df[u'Y坐标'] = df[u'Y坐标'].map(lambda x: x.replace('Y:', ''))

    #度分秒转换为经纬度坐标
    df[u'X坐标'] = df[u'X坐标'].apply(longitude_, )
    df[u'Y坐标'] = df[u'Y坐标'].apply(longitude_, )

    df.to_excel(u'C:\\Users\\YaoJun\\Desktop\\矿产图斑列表1.xlsx', index=False)

    #创建输出shp文件
    shpfname = xlsx_path.split('.')[0] + ".shp"
    CreatePolygonShp(shpfname)

    #添加所有列名为字段
    for col_name in df.columns[:-2]:
        if len(col_name) > 5:
            continue
        if u"面积" in col_name:
            AddField(shpfname, col_name, ogr.OFTReal)
        else:
            AddField(shpfname, col_name, ogr.OFTString, 30)

    #打开空白shp
    driver = ogr.GetDriverByName('ESRI Shapefile')
    fc1 = driver.Open(shpfname, 1)
    layer = fc1.GetLayer()
    layer_defn = layer.GetLayerDefn()
    field_count = layer_defn.GetFieldCount()

    #添加图形与属性信息
    for i in range(0, len(df.index)):
        feature = ogr.Feature(layer.GetLayerDefn())
        #依次给字段附属性值
        for j in range(field_count):
            field_defn = layer_defn.GetFieldDefn(j)
            field_name = field_defn.GetName()
            if field_name in df.columns:
                feature.SetField(field_name, str(df.iloc[i][field_name]))
        #创建点
        wkt = "POINT({} {})".format(df.iloc[i][u'X坐标'], df.iloc[i][u'Y坐标'])  # "point({} {})"两个参数之间没有逗号，记住
        point = ogr.CreateGeometryFromWkt(wkt)
        #print(point)
        feature.SetGeometry(point)
        layer.CreateFeature(feature)
        del feature

    pass

4.实验结果