【GIS】GDAL Python 影像裁剪

# -*- coding: utf-8 -*-
"""
Created on Fri Nov 30 11:45:03 2018

@author: Administrator
"""

from osgeo import gdal
from osgeo import osr
import numpy as np
import math
import time

lonMeter = 0.00001141 
latMeter = 0.00000899

#MeterParam = 0.00001 * 42496 / (124.44282531738276-124.3288421630859)
MeterParam = 3.7282702222226876

def getSRSPair(dataset):
    '''
    獲得給定數據的投影參考系和地理參考系
    :param dataset: GDAL地理數據
    :return: 投影參考系和地理參考系
    '''
    prosrs = osr.SpatialReference()
    prosrs.ImportFromWkt(dataset.GetProjection())
    geosrs = prosrs.CloneGeogCS()
    return prosrs, geosrs

def geo2lonlat(dataset, x, y):
    '''
    將投影坐標轉為經緯度坐標（具體的投影坐標系由給定數據確定）
    :param dataset: GDAL地理數據
    :param x: 投影坐標x
    :param y: 投影坐標y
    :return: 投影坐標(x, y)對應的經緯度坐標(lon, lat)
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(prosrs, geosrs)
    coords = ct.TransformPoint(x, y)
    return coords[:2]


def lonlat2geo(dataset, lon, lat):
    '''
    將經緯度坐標轉為投影坐標（具體的投影坐標系由給定數據確定）
    :param dataset: GDAL地理數據
    :param lon: 地理坐標lon經度
    :param lat: 地理坐標lat緯度
    :return: 經緯度坐標(lon, lat)對應的投影坐標
    '''
    prosrs, geosrs = getSRSPair(dataset)
    ct = osr.CoordinateTransformation(geosrs, prosrs)
    coords = ct.TransformPoint(lon, lat)
    return coords[:2]

def imagexy2geo(dataset, row, col):
    '''
    根據GDAL的六參數模型將影像圖上坐標（行列號）轉為投影坐標或地理坐標（根據具體數據的坐標系統轉換）
    :param dataset: GDAL地理數據
    :param row: 像素的行號
    :param col: 像素的列號
    :return: 行列號(row, col)對應的投影坐標或地理坐標(x, y)
    '''
    trans = dataset.GetGeoTransform()
    px = trans[0] + col * trans[1] + row * trans[2]
    py = trans[3] + col * trans[4] + row * trans[5]
    return px, py


def geo2imagexy(dataset, x, y):
    '''
    根據GDAL的六 參數模型將給定的投影或地理坐標轉為影像圖上坐標（行列號）
    :param dataset: GDAL地理數據
    :param x: 投影或地理坐標x
    :param y: 投影或地理坐標y
    :return: 影坐標或地理坐標(x, y)對應的影像圖上行列號(row, col)
    '''
    trans = dataset.GetGeoTransform()
    a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
    b = np.array([x - trans[0], y - trans[3]])
    return np.linalg.solve(a, b)  # 使用numpy的linalg.solve進行二元一次方程的求解

def imagexy2lonlat(dataset,row, col):
    '''
    影像行列轉經緯度：
    ：通過影像行列轉平面坐標
    ：平面坐標轉經緯度
    '''
    coords = imagexy2geo(dataset, row, col)
    coords2 = geo2lonlat(dataset,coords[0], coords[1])
    return (coords2[0], coords2[1])

def lonlat2imagexy(dataset,x, y):
    '''
    影像行列轉經緯度：
    ：通過經緯度轉平面坐標
    ：平面坐標轉影像行列
    '''
    coords = lonlat2geo(dataset, x, y)
    coords2 = geo2imagexy(dataset,coords[0], coords[1])
    return (int(round(abs(coords2[0]))), int(round(abs(coords2[1]))))


if __name__ == '__main__':
    gdal.AllRegister()
    dataset = gdal.Open(r"D:\RSData\DAQING_SHAERTU\薩爾圖區_大圖：拼接\L19.tif")
    
    print('數據投影：')
    projection = dataset.GetProjection()
    print(projection)
    
    print('數據的大小（行，列）：')
    print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))
    
    geotransform = dataset.GetGeoTransform()
    print('地理坐標：')
    print(geotransform)
 
    x = 464201
    y = 5818760
    lon = 122.47242
    lat = 52.51778
    row = 0
    col = 0

#    print('投影坐標 -> 經緯度：')
#    coords = geo2lonlat(dataset, x, y)
#    print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
#    
#    print('經緯度 -> 投影坐標：')
#    coords = lonlat2geo(dataset, lon, lat)
#    print('(%s, %s)->(%s, %s)' % (lon, lat, coords[0], coords[1]))
#
#    print('圖上坐標 -> 投影坐標：')
#    coords = imagexy2geo(dataset, row, col)
#    print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1]))
#    
#    print('投影坐標 -> 圖上坐標：')
#    coords = geo2imagexy(dataset, x, y)
#    print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
    
#    print('圖上坐標 -> 投影坐標：')
#    coords = imagexy2geo(dataset, row, col)
#    print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1]))
#    print('投影坐標 -> 經緯度：')
#    coords2 = geo2lonlat(dataset,coords[0], coords[1])
#    print('(%s, %s)->(%s, %s)' % (coords[0], coords[1], coords2[0], coords2[1]))
    
#    coords = imagexy2lonlat(dataset, row, col)
#    print('影像像素 -> 經緯度：')
#    print('(%s, %s)->(%s, %s)' % ( row, col, coords[0], coords[1]))
#    coords = imagexy2lonlat(dataset, dataset.RasterXSize, dataset.RasterYSize)
#    print('影像像素 -> 經緯度：')
#    print('(%s, %s)->(%s, %s)' % ( dataset.RasterXSize, dataset.RasterYSize, coords[0], coords[1]))
#    
#    coords = lonlat2imagexy(dataset, 124.3288421630859, 46.391464001559044)
#    print('經緯度 -> 影像像素 ：')
#    print('(%s, %s)->(%s, %s)' % ( 124.3288421630859, 46.391464001559044, coords[0], coords[1]))
#    coords = lonlat2imagexy(dataset, 124.44282531738276, 46.32796494040744)
#    print('經緯度 -> 影像像素 ：')
#    print('(%s, %s)->(%s, %s)' % ( 124.44282531738276, 46.32796494040744, coords[0], coords[1])) 
    
    #經緯度轉像素
    xoffset=0
    yoffset=0
    
    x,y = 125.059,46.894
    
    xoffset,yoffset = lonlat2imagexy(dataset, x,y)
    print('坐標轉換-對應行列像素位置')
    print('(%s, %s)->(%s, %s)' % (x,y, xoffset,yoffset))
    
    width=int(500 * MeterParam)
    height=int(500 * MeterParam)
    
    
    if xoffset - width <= 0  and yoffset - height <= 0 :
        print("左上角")
        xoffset = 0
        yoffset = 0
    elif xoffset - width <= 0  and yoffset - height > 0 :
        print("左邊")
        xoffset = 0
    elif xoffset - width > 0  and yoffset - height <= 0 :
        print("頂邊")
        yoffset = 0      
    else :
        print("中間區域")
        xoffset = xoffset - width
        yoffset = yoffset - height

    
    width = width * 2
    height = height * 2
    
    print('切割范圍')
    print('寬高(%s, %s)->偏移起點(%s, %s)' % (width, height, xoffset,yoffset))
#    xoffset,yoffset,width,height = 175360/2,123136/2,1000,1000
    
    newData = np.zeros([width,height,3])
    band = dataset.GetRasterBand(1) 
    r=band.ReadAsArray(xoffset,yoffset,width,height)
    NoData = band.GetNoDataValue()
    newData[:,:,0] = r
    
    band = dataset.GetRasterBand(2)
    g=band.ReadAsArray(xoffset,yoffset,width,height)
    
    band = dataset.GetRasterBand(3)
    b=band.ReadAsArray(xoffset,yoffset,width,height)
    
    ticks = time.time()
    resultPath = "D:\\RS%s.jpg" % (int(ticks))
    
    newData[:,:,0] = r
    newData[:,:,1] = g
    newData[:,:,2] = b
    
    format = "GTiff"   
    driver = gdal.GetDriverByName(format)
    ds = driver.Create(resultPath, width, height, 3, gdal.GDT_Float32)
    geotransform1 = geotransform
    px = geotransform[0] + xoffset * geotransform[1] + yoffset * geotransform[2]
    py = geotransform[3] + xoffset * geotransform[4] + yoffset * geotransform[5]
    geotransform1 = (px, 0.29858214173896974, 0.0, py, 0.0, -0.29858214173896974)
#    print(geotransform1[0])
    ds.SetGeoTransform(geotransform1)
    ds.SetProjection(projection)
    lay01= ds.GetRasterBand(1)
    lay02= ds.GetRasterBand(2)
    lay03= ds.GetRasterBand(3)
#    ds.GetRasterBand(1).SetNoDataValue(0)
#    ds.GetRasterBand(2).SetNoDataValue(0)
#    ds.GetRasterBand(3).SetNoDataValue(0)
    lay01.WriteArray(b) 
    lay02.WriteArray(g)
    lay03.WriteArray(r)
#    ds.FlushCache()
#    ds = None
    del ds
 
    
    import cv2

    import matplotlib.pyplot as plt
    
    
    img2=cv2.merge([r,g,b])
    plt.imshow(img2)
    plt.xticks([]),plt.yticks([]) # 不顯示坐標軸
    plt.show()
    
    ticks = time.time()
#    cv2.imwrite("D:\\RS%s.jpg" % (int(ticks)) , img2)
    
    print("OK")
 
    

 轉自：https://blog.csdn.net/theonegis/article/details/54427906