当前位置: 代码迷 >> 综合 >> 【GDAL】python批量读取文件夹下.gz压缩文件并根据压缩文件名称生成点状shp
  详细解决方案

【GDAL】python批量读取文件夹下.gz压缩文件并根据压缩文件名称生成点状shp

热度:68   发布时间:2023-11-30 17:04:11.0

      之前写过一篇【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 filenamefilename = []  # 文件名pointListX = []  # X坐标pointListY = []  # Y坐标Date = []  # 日期ProductID = []  # 产品IDSensorType = []  # 传感器类型# 设置中文字符编码gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")gdal.SetConfigOption("SHAPE_ENCODING", "GBK")# register a driverdriver = ogr.GetDriverByName('ESRI Shapefile')# outfile nameoutfile = '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 fileoutDS = driver.CreateDataSource(outfile)if outDS is None:print('Cannot open this file: ', outfile)# define projectionproj = osr.SpatialReference()proj.ImportFromEPSG(4326)# create layeroutLayer = outDS.CreateLayer('point', proj, geom_type=ogr.wkbPoint)# create fieldsfieldDefn1 = 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 defintionoutFeatureDefn = outLayer.GetLayerDefn()# write every point to shape filefor 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!')