python之tiff

1、安装libtiff

conda安装的没生效

libtiff                   4.1.0                h56a325e_0    defaults

下载whl https://www.lfd.uci.edu/~gohlke/pythonlibs/#pylibtiff

pip install pylibtiff-0.4.2-cp37-cp37m-win_amd64.whl

pylibtiff                 0.4.2                    pypi_0    pypi

2、读取tiff数据

from libtiff import TIFF

    tiff_image_name = r'..\tif\baoyu\baoyu.tif'
    tif = tiff2Stack(tiff_image_name)
    print(tif)
    
def tiff2Stack(filePath):
    tif = TIFF.open(filePath,mode='r')
    stack = []
    for img in list(tif.iter_images()):
        stack.append(img)
    return stack

conda install tifffile

import tifffile

tif = tifffile.imread(tiff_image_name)

def tiff_stack(tiff_image_name, file_time_str, type):
    with tifffile.TiffFile(tiff_image_name) as tif:
        print(tif)
        meta_data = tif.geotiff_metadata
        xDelta = meta_data['ModelPixelScale'][0]
        yDelta = meta_data['ModelPixelScale'][1]
        start_lon = meta_data['ModelTiepoint'][3]
        end_lat = meta_data['ModelTiepoint'][4]
        data = tif.asarray()
        ImageLength = data.shape[0]
        ImageWidth = data.shape[1]
        end_lon = start_lon + ImageWidth * xDelta
        start_lat = end_lat - ImageLength * yDelta
        dat = data.reshape(1, 1, ImageLength, ImageWidth)
        dat[dat == -9999] = np.nan
        lons, lats = [], []
        for x in range(ImageWidth):
            lons.append(start_lon + x * xDelta)
        for y in range(ImageLength):
            lats.append(end_lat - y * yDelta)
       	save_file = Path(tiff_image_name).parent.joinpath(f'{type}_{file_time_str}.zip')
        data_to_upload(dat, lats, lons, file_time_str, save_file, 'CQIC_GRID_TIFF', type)

3、shp截取tif

1.shp必须只有1个范围

from osgeo import gdal

    input_raster = r"E:/chqqh/tif/FA/o0001.tif"
    input_raster = gdal.Open(input_raster)
    input_shape = r"E:/chqqh/tif/shp/1.shp"  # or any other format
    output_raster = r'E:\chqqh\tif\test\o0001_1.tif'
    ds = gdal.Warp(output_raster,
                   input_raster,
                   format='GTiff',
                   cutlineDSName=input_shape,  # or any other file format
                   cutlineWhere="FIELD = 'whatever'",
                   # optionally you can filter your cutline (shapefile) based on attribute values
                   dstNodata=-9999)  # select the no data value you like
    ds = None

4、获取tif的经纬范围

from osgeo import gdal
import glob

min_x, max_y, max_x, min_y = get_extent(f'E:/tif/test/1.tif')

def get_extent(in_fn):
    ds = gdal.Open(in_fn)
    geotrans = list(ds.GetGeoTransform())
    xsize = ds.RasterXSize
    ysize = ds.RasterYSize
    min_x = geotrans[0]
    max_y = geotrans[3]
    max_x = geotrans[0] + xsize * geotrans[1]
    min_y = geotrans[3] + ysize * geotrans[5]
    ds = None
    return min_x, max_y, max_x, min_y

5、多个tif合并

from osgeo import gdal
import glob
import rasterio as rio
from rasterio.merge import merge
from rasterio.plot import show
 
    in_files = glob.glob(r'E:/chqqh/tif/test/*.tif')
    src_files_to_merged = []
    for fp in in_files:
        tif = rio.open(fp)
        src_files_to_merged.append(tif)
    stiched, out_trans = merge(src_files_to_merged)
    # show(stiched, cmap='terrain')
    out_meta = tif.meta.copy()
    out_meta.update({"driver": "GTiff", "height": stiched.shape[1], "width": stiched.shape[2], "transform": out_trans})
    outfn = r'E:/chqqh/tif/fa_result/o0001_res.tif'
    with rio.open(outfn, "w", **out_meta) as dest:
        dest.write(stiched)

6、nparray生成tiff

import osgeo.osr as osr
from osgeo import gdal
import numpy as np

    rect_json = net_io.rect_post('106.54', '30.16', '105.8', '29.59', 'TMP', '2020-11-12', '08', '24')
    res = rect_json['gridDatas'][0]
    values = res['values'] # type(list)
    startLat = res['startLat'] #
    startLon = res['startLon']
    endLat = res['endLat']
    endLon = res['endLon']
    cols = res['cols'] #30
    rows = res['rows'] #23
    element = res['elements']
    valid = res['valid']
    latResolution = float(res['latResolution'])
    lonResolution = float(res['lonResolution'])
    dat = np.array(values).astype("float32") #.reshape(1, 1, rows, cols)
    # 创建.tif文件
    driver = gdal.GetDriverByName('GTiff')
    out_tif_name = r'E:\test\TMP.tif'
    out_tif = driver.Create(out_tif_name, cols, rows, 1, gdal.GDT_Float32)  # 创建框架

    # 设置影像的显示范围
    # Lat_Res一定要是-的
    geotransform = (startLon, lonResolution, 0, endLat, 0, -latResolution)
    out_tif.SetGeoTransform(geotransform)

    # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
    out_tif.SetProjection(srs.ExportToWkt())  # 给新建图层赋予投影信息

    # 数据写出
    out_tif.GetRasterBand(1).WriteArray(dat[::-1])  # 将数据写入内存，此时没有写入硬盘 数据维度是正序排列，需要逆序一下
    out_tif.FlushCache()  # 将数据写入硬盘
    out_tif = None  # 注意必须关闭tif文件