python gdal多值提取至点工具

提取shp点文件栅格目录中的栅格数值序列

import os
from osgeo import gdal, gdalconst, ogr
import pandas as pd

class ExtractByPoint(object):
    def __init__(self, x, y, tifFile) -> None:
        self.x = x
        self.y = y
        self.tifFile = tifFile

# 地理坐标转像素坐标
    @staticmethod
    def geo2pixel(geot: list, geo: list):
        px = (geot[5]*geo[0]-geot[2]*geo[1]-geot[0]*geot[5] +
              geot[2]*geot[3])/(geot[1]*geot[5]-geot[2]*geot[4])
        py = (geot[4]*geo[0]-geot[1]*geo[1]+geot[1]*geot[3] -
              geot[4]*geot[0])/(geot[2]*geot[4]-geot[1]*geot[5])
        return [int(px), int(py)]

    def getValue(self):
        ds = gdal.Open(self.tifFile, gdalconst.GA_ReadOnly)
        geot = gdal.Dataset.GetGeoTransform(ds)
        px, py = ExtractByPoint.geo2pixel(geot, [self.x, self.y])
        band = gdal.Dataset.GetRasterBand(ds, 1)
        bandArray = gdal.Band.ReadAsArray(band)
        value = bandArray[py, px]
        return format(value, "0.2f")


if __name__ == "__main__":
    # 提取所有的tif影像
    tifsdir = 'tifs'
    alltifs = dict()
    for tifsdir_ in os.listdir(tifsdir):
        tifs = os.listdir(os.path.join(tifsdir, tifsdir_))
        tifs=filter(lambda filename:filename.endswith('.tif'),tifs)
        tifs = [os.path.join(tifsdir, tifsdir_, tif) for tif in tifs]
        alltifs[tifsdir_] = tifs
    pointShp = 'cape\cape.shp'
    ogr.RegisterAll()
    pds = ogr.Open(pointShp)
    pLay = ogr.DataSource.GetLayerByIndex(pds, 0)
    ogr.Layer.ResetReading(pLay)
    # 便利所有的点
    while True:
        pFea = ogr.Layer.GetNextFeature(pLay)
        if not(pFea):
            break
        name = int(ogr.Feature.GetFieldAsDouble(pFea, "ForecastSt"))
        pGeo = ogr.Feature.GetGeometryRef(pFea)
        x = ogr.Geometry.GetX(pGeo)
        y = ogr.Geometry.GetY(pGeo)
        pdict = dict()
        # 添加日期字段
        for key in alltifs.keys():
            tifs=alltifs[key]
            for tif in tifs:
                tifname=os.path.basename(tif)
                index=tifname.find('_')
                dates.append(tifname[(index+1):-4].replace('_',''))
            break
        for key, tifs in alltifs.items():
            pvalues = []
            for tif in tifs:
                e = ExtractByPoint(x, y, tif)
                pvalue = e.getValue()
                pvalues.append(pvalue)
            pdict[key] = pvalues
        dates = []
        pdict['date']=dates
        pdf = pd.DataFrame(pdict)
        csv_name = './csvs/{}.csv'.format(name)
        pdf.to_csv(csv_name, index=False)
    
    print('ok!')