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