个人整理的python地理信息库,本文代码也在其中，欢迎使用: https://github.com/Adam0429/geo-py
顺便帮公司打个广告，欢迎遥感类专家的加入～～
投简历戳这里，英视睿达欢迎您！

可能安装依赖时会遇到问题，本人提供有偿远程安装服务，需要的可以私信

废话不多说，直接上代码

from tqdm import tqdm
from osgeo import gdal,ogr,osr
import numpy as np
from glob import glob

def read_tif(tif_path):
    ds = gdal.Open(tiff_path)
    row = ds.RasterXSize
    col = ds.RasterYSize
    band = ds.RasterCount

    for i in range(band):
        data = ds.GetRasterBand(i+1).ReadAsArray()

        data = np.expand_dims(data , 2)
        if i == 0:
            allarrays = data
        else:
            allarrays = np.concatenate((allarrays, data), axis=2)
    return {'data':allarrays,'transform':ds.GetGeoTransform(),'projection':ds.GetProjection(),'bands':band,'width':row,'height':col}
    # 左上角点坐标 GeoTransform[0],GeoTransform[3] Transform[1] is the pixel width, and Transform[5] is the pixel height

def write_tif(fn_out, im_data, transform,proj=None):
    '''
    功能:
    ----------
    将矩阵按某种投影写入tif，需指定仿射变换矩阵，可选渲染为rgba
    
    参数:
    ----------
    fn_out:str
        输出tif图的绝对文件路径
    im_data: np.array
        tif图对应的矩阵
    transform: list/tuple 
        gdal-like仿射变换矩阵，若im_data矩阵起始点为左上角且投影为4326，则为
            (lon_x.min(), delta_x, 0, 
             lat_y.max(), 0, delta_y)
    proj: str（wkt格式）
        投影，默认投影坐标为4326，可用osr包将epsg转化为wkt格式，如
            srs = osr.SpatialReference()# establish encoding
            srs.ImportFromEPSG(4326)    # WGS84 lat/lon
            proj = srs.ExportToWkt()    # create wkt fromat of proj


    # 设置投影，proj为wkt format
    if proj is None:
        proj = 'GEOGCS["WGS 84",\
                     DATUM["WGS_1984",\
                             SPHEROID["WGS 84",6378137,298.257223563, \
                                    AUTHORITY["EPSG","7030"]], \
                             AUTHORITY["EPSG","6326"]], \
                     PRIMEM["Greenwich",0, \
                            AUTHORITY["EPSG","8901"]], \
                     UNIT["degree",0.0174532925199433, \
                            AUTHORITY["EPSG","9122"]],\
                     AUTHORITY["EPSG","4326"]]'
    # 渲染为rgba矩阵
    # 设置数据类型
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
    # 将(通道数、高、宽)顺序调整为(高、宽、通道数)
    # print('shape of im data:', im_data.shape)
    im_bands = min(im_data.shape)
    im_shape = list(im_data.shape)
    im_shape.remove(im_bands)
    im_height, im_width = im_shape
    band_idx = im_data.shape.index(im_bands)
    # 找出波段是在第几个

    # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(fn_out, im_width, im_height, im_bands, datatype)

    # if dataset is not None:
    dataset.SetGeoTransform(transform)  # 写入仿射变换参数
    dataset.SetProjection(proj)  # 写入投影

    if im_bands == 1:

        # print(im_data[:, 0,:].shape)
        if band_idx == 0:
            dataset.GetRasterBand(1).WriteArray(im_data[0, :, :])
        elif band_idx == 1:
            dataset.GetRasterBand(1).WriteArray(im_data[:, 0, :])
        elif band_idx == 2:
            dataset.GetRasterBand(1).WriteArray(im_data[:, :, 0])

    else:

        for i in range(im_bands):
            if band_idx == 0:
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i, :, :])
            elif band_idx == 1:
                dataset.GetRasterBand(i + 1).WriteArray(im_data[:, i, :])
            elif band_idx == 2:
                dataset.GetRasterBand(i + 1).WriteArray(im_data[:, :, i])

    dataset.FlushCache()
    del dataset
    driver = None

使用案例

>>> geo.read_tif('/Users/wangfeihong/Desktop/GF2_PMS2_E102.4_N26.3_20210420_L1A0005608668-MSS2_clip.tif')
{'data': array([[[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        ...,
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]],

       [[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        ...,
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]],

       [[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        ...,
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]],

       ...,

       [[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        ...,
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]],

       [[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        ...,
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]],

       [[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        ...,
        [0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]]], dtype=uint16), 'transform': (102.36514511007022, 7.796940285642904e-06, 0.0, 26.32195006415771, 0.0, -7.797041881219687e-06), 'projection': 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]', 'bands': 4, 'width': 11577, 'height': 6836}