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

      之前写过一篇【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!')

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值