1. GetGeoTransform()使用、gdal截取图像

GetGeoTransform()
GeoTransform[0]，左上角横坐标（应该是投影坐标）
GeoTransform[2]，行旋转
GeoTransform[1]，像元宽度(影像在水平空间的分辨率)
GeoTransform[3]，左上角纵坐标（应该是投影坐标）
GeoTransform[4]，列旋转
GeoTransform[5]，像元高度(影像在垂直空间的分辨率)，
如果影像是指北的，GeoTransform[2]和GeoTransform[4]这两个参数的值为0，GeoTransform[5]为负；
如果图像不含地理坐标信息，默认返回值是：(0,1,0,0,0,1)

ReadAsArray(xoff,yoff,xsize,ysize)
#xoff,yoff是取值窗口的左上角在实际数据中所处象元的xy位置。
#xsize,ysize是取值窗口覆盖的区域大小

import gdal  #导入库

dataset = gdal.open("filename")  #文件名，如*.tif
geotransform = dataset.GetGeoTransform()
print('数据投影：')
print(dataset.GetProjection())
width = int(dataset.RasterXSize) #读取图像的宽度，x方向上的像素个数
height = int(dataset.RasterYSize)  #读取图像的高度，y方向上的像素个数
RasterCount = int(dataset.RasterCount)
print('数据的大小（行/height，列/weight, RasterCount）：')
print('(%s %s %s)' % (height, width, RasterCount))

if geotransform[0] == 0.0 and geotransform[1] == 1.0:
	print("EMPTY IMAGE!")
img_tu = dataset.ReadAsArray(xoff=xp_1, yoff=yp_1, xsize=xp_2 - xp_1, ysize=yp_2 - yp_1) #截取原图像中左上角(xp_1,yp_1)至右下角(xp_1,yp_1)之间的矩形区域

# 之后就可对截取区域进行操作/直接存储到路径了

2. 使用GDAL进行影像投影坐标、地理坐标、图上坐标的转换

转载来源：https://theonegis.blog.csdn.net/article/details/54427906?spm=1001.2014.3001.5506
亲测好用！码一下，侵删。

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

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(lat, lon)
    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进行二元一次方程的求解


if __name__ == '__main__':
    gdal.AllRegister()
    dataset = gdal.Open(r"F:\2016\Data\Great Khingan\DEM\Projection\strm_6102_UTM.tif")
    print('数据投影：')
    print(dataset.GetProjection())
    print('数据的大小（行，列）：')
    print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))

    x = 464201
    y = 5818760
    lon = 122.47242
    lat = 52.51778
    row = 2399
    col = 3751

    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]))