之前写过一篇【GDAL】python读取txt影像名称文件生成shp,本篇博客是对之前那篇的扩展。之前那篇是把所有.gz文件名写在了一个txt文件里,通过读取txt文件将里面的文件名保存到一个文件名列表中即可。但是这篇博客是读取一个文件夹下所有文件夹中的所有.gz压缩包,将这些压缩包的文件路径名存放在一个列表中解析出文件名路径、XY、传感器类型、日期和产品ID,创建一个点状shape文件,并将这些属性保存到属性表中
用到的主要是GDAL库
文件夹下文件如下图样子:
代码运行结果如下图所示
代码如下:
import gdal
from gdal import ogr
from gdal import osr
import os
# read txt and return a filename list
def readFile(filepath):
# 若filepath是一个文件路径,则直接读取该文件
# f = open(filename, 'r')
# lines = f.readlines()
# pointFilenameList = []
#
# for line in lines:
# pointFilenameList.append(line)
#
# return pointFilenameList
# 获得该路径下的所有文件夹
list = os.listdir(filepath)
# 读所有文件夹下所有gz文件
for item in list:
tmp_path = os.path.join(filepath, item)
if not os.path.isdir(tmp_path):
if item.split('.')[-1] == 'gz':
gzlist.append(tmp_path.split(':')[1])
else:
readFile(tmp_path)
# 读单个文件夹
# for item in list:
# if item.split('.')[-1] == 'gz':
# gzlist.append(item)
return gzlist
# input filename list and create point shape file
def createShp(pointFilenameList):
# define list to restore x y coor and filename
filename = [] # 文件名
pointListX = [] # X坐标
pointListY = [] # Y坐标
Date = [] # 日期
ProductID = [] # 产品ID
SensorType = [] # 传感器类型
# 设置中文字符编码
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
# register a driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# outfile name
outfile = 'GF.shp'
for line in pointFilenameList:
print(line)
nameline = line.split('\\')[-1]
sensor = nameline.split('.')[-3].split('-')[-1][0:-1]
if nameline.split('_')[0][0:-1] == 'GF': # GF数据
if sensor == 'MSS':
y = nameline.split('_')[2][1:]
x = nameline.split('_')[3][1:]
date = nameline.split('_')[4]
pid = str(int(nameline.split('_')[5].split('-')[0].split('A')[1]))
stype = nameline.split('_')[0]
pointListX.append(x)
pointListY.append(y)
SensorType.append(stype)
Date.append(date)
ProductID.append(pid)
filename.append(line)
elif line.split('_')[0][0:-1] == 'ZY': # ZY数据
y = nameline.split('_')[2][1:]
x = nameline.split('_')[3][1:]
date = nameline.split('_')[4]
pid = str(int(nameline.split('_')[5].split('.')[0].split('A')[1]))
stype = nameline.split('_')[0]
pointListX.append(x)
pointListY.append(y)
SensorType.append(stype)
Date.append(date)
ProductID.append(pid)
filename.append(line)
# # GF3
# y = line.split('_')[4][1:]
# x = line.split('_')[5][1:]
# date = line.split('_')[6]
# pid = str(int(line.split('_')[9].split('.')[0][1:]))
# stype = line.split('_')[0]
#
# pointListX.append(x)
# pointListY.append(y)
# SensorType.append(stype)
# Date.append(date)
# ProductID.append(pid)
# filename.append(line)
if os.path.exists(outfile):
driver.DeleteDataSource(outfile)
# create shape file
outDS = driver.CreateDataSource(outfile)
if outDS is None:
print('Cannot open this file: ', outfile)
# define projection
proj = osr.SpatialReference()
proj.ImportFromEPSG(4326)
# create layer
outLayer = outDS.CreateLayer('point', proj, geom_type=ogr.wkbPoint)
# create fields
fieldDefn1 = ogr.FieldDefn('filename', ogr.OFTString)
fieldDefn1.SetWidth(60)
outLayer.CreateField(fieldDefn1, 1)
fieldDefn2 = ogr.FieldDefn('X', ogr.OFTString)
fieldDefn2.SetWidth(4)
outLayer.CreateField(fieldDefn2, 1)
fieldDefn3 = ogr.FieldDefn('Y', ogr.OFTString)
fieldDefn3.SetWidth(5)
outLayer.CreateField(fieldDefn3, 1)
fieldDefn4 = ogr.FieldDefn('Sensor', ogr.OFTString)
fieldDefn4.SetWidth(5)
outLayer.CreateField(fieldDefn4, 1)
fieldDefn5 = ogr.FieldDefn('Date', ogr.OFTString)
fieldDefn5.SetWidth(8)
outLayer.CreateField(fieldDefn5, 1)
fieldDefn6 = ogr.FieldDefn('ProductID', ogr.OFTString)
fieldDefn6.SetWidth(10)
outLayer.CreateField(fieldDefn6, 1)
# get feature defintion
outFeatureDefn = outLayer.GetLayerDefn()
# write every point to shape file
for i in range(len(pointListX)):
outFeature = ogr.Feature(outFeatureDefn)
wkt = "POINT(%f %f)" % (float(pointListY[i]), float(pointListX[i]))
point = ogr.CreateGeometryFromWkt(wkt)
outFeature.SetGeometry(point)
outFeature.SetField('filename', filename[i])
outFeature.SetField('X', pointListX[i])
outFeature.SetField('Y', pointListY[i])
outFeature.SetField('Sensor', SensorType[i])
outFeature.SetField('Date', Date[i])
outFeature.SetField('ProductID', ProductID[i])
outLayer.CreateFeature(outFeature)
outFeature.Destroy()
outDS.Destroy()
if __name__ == '__main__':
filepath = 'G:/GF'
gzlist = []
os.chdir(filepath)
namelist = readFile(filepath)
createShp(namelist)
print('Success!')