python/gdal处理遥感影像（读取、投影转换、裁剪、建立图像金字塔等）

gdal库简单介绍python使用gdal一、安装python环境二、安装gdal库三、使用gdal处理遥感影像1、读取TIFF数据2、获取波段信息和波段数据3 写出tif文件4、金字塔生成5 投影转换

gdal库简单介绍

gdal的全称含义是地理空间数据抽象库，是一个在X/MIT许可协议下读写空间数据（栅格数据和矢量数）的开源库，源码编写使用c/c++，但也可以在 Perl、python、VB6、Ruby、Java、C# 等语言中调用 GDAL,这令 GDAL 的应用变得非常广泛。目前国内广泛用于遥感图像的处理平台GIS和RS的底层都是使用GDAL来读写空间数据。

python使用gdal

一、安装python环境

请移步百度搜索，python的安装。

二、安装gdal库

gdal库依赖于多个库，为了避免不同库之间的版本冲突，新建一个虚拟环境，来安装gdal库最好。

使用conda安装

conda install gdal -c conda-forge

使用pip安装

使用pip安装，可参照上面conda的方式，也可以在python扩展包网站上下载对应python环境版本的whl文件进行安装。python扩展网站:Unofficial Windows Binaries for Python Extension Packages

pip install c:/download/GDAL-3.2.3-cp37-cp37m-win_amd64.whl

搜索结果如下：

三、使用gdal处理遥感影像

1、读取TIFF数据

import gdal import osrfrom osgeo import gdalgdal.AllRegister() #先载入数据驱动，也就是初始化一个对象，让它“知道”某种数据结构，但是只能读，不能写ds = gdal.Open("d:/test/test.tif") # 打开文件bands = ds.RasterCount()# 获取波段数img_width,img_height = ds.RasterXSize,ds.RasterYSize # 获取影像的宽高geotrans = ds.GetGeoTransform() # 获取影像的投影信息

投影信息说明

在GDAL中，这六个参数包括左上角坐标（即图像的经度（lon），维度（lat)，像元X、Y方向大小，旋转等信息。 要注意，Y方向的像元大小为负值。六参数的顺序为（lon,res,0,lat,-res,0）

2、获取波段信息和波段数据

band = ds.GetRasterBand(1) #获取波段信息,如有多个波段，数字递增，每一个数字表示第几个波段的信息。xsize.ysize = band.XSize,band.Ysizedtye = band.DataTypenodata = band.GetNoDataValue()maxdata = band.GetMaximum() # 不包括nodatabdata = band.ReadAsArray()

ReadAsArray(self, xoff=0, yoff=0, xsize=None, ysize=None, buf_obj=None, buf_xsize = None, buf_ysize = None, buf_type = None,resample_alg = GRIORA_NearestNeighbour, callback = None, callback_data = None):

看看函数的几个参数的意义。头两个100是取值窗口的左上角在实际数据中所处象元的xy位置。后两个是取值窗口覆盖的区域大小，再后面 是取值窗口取出数组进行缩放后数组的大小。这里需要注意的是这里的buffer大小是根据参数自动分配的，可以不指定，如果不指定，则和第3，4个参数一致。经过5，6两个参数的设置，可以进行缩放。resample_alg 则表示缩放的方式，即重采样方式。

3 写出tif文件

import netCDF4 as ncimport gdal,osrncfile = NCOper('F:\FY3D_MERSI_ORBD_L2_NCC_MLT_GLL_2001_POAD_5000M_V0.HDF')output = 'F:\FY3D_MERSI_ORBD_L2_NCC_MLT_GLL_2001_POAD_5000M_V0.HDF.tif'dat = nc.Dataset(ncfile).variables['NCC'][:]col,row = dat.shapedriver = gdal.GetDriverByName('GTiff')dst_ds = driver.Create(output, row, col,1,gdal.GDT_Float32)dst_ds.SetGeoTransform([-180,0.05,0.0,90,0.0,-0.05])dst_ds.GetRasterBand(1).WriteArray(dat)proj = osr.SpatialReference()proj.SetWellKnownGeogCS('WGS84')# 定义参考dst_ds.SetProjection(proj.ExportToWkt())dst_ds.FlushCache()

4、金字塔生成

图像金子塔是一种多分辨率图像组合的结构。底部是待处理图像的高分辨率表示，顶部是低分辨率的近似。当从下层向上移动时，尺寸和分辨率都会下降。list里面是要生成金字塔的层级，全部为2的幂次，第一层是原始数据。利用gdal生成的ovr数据，实际上就是一个多分辨率的tif文件。文件明在原始的tif后面直接加追.ovr。

ds.BuildOverviews(overviewlist=[1,2,4,8,16])

5 投影转换

gdal投影转换参考文章链接: 投影转换，

下面这个投影是基于pyproj进行的

from pyproj import Proj, transform,Transformer# 方式一ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')x, y, z = pyproj.transform(lla, ecef, lon, lat, height, radians=False) # radians否用弧度返回值# 方式一随着proj版本的升级，可能会被弃用，推荐方式二# 方式二transprojr = Transformer.from_crs("EPSG:4326",{"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},always_xy=True,)xpjr, ypjr, zpjr = transprojr.transform(lon,lat,height, radians=False)

以下内容转自 /article/28411277112/

# -*- encoding: utf-8 -*-from osgeo import gdalfrom osgeo import osrimport numpy as npdef getSRSPair(dataset):'''获得给定数据的投影参考系和地理参考系:param dataset: GDAL地理数据:return: 投影参考系和地理参考系'''prosrs = osr.SpatialReference()prosrs.ImportFromWkt(dataset.GetProjection())geosrs = prosrs.CloneGeogCS()return prosrs, geosrsdef 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, pydef 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:\\Data\Great Khingan\DEM\Projection\strm_6102_UTM.tif")print('数据投影：')print(dataset.GetProjection())print('数据的大小（行，列）：')print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))x = 464201y = 5818760lon = 122.47242lat = 52.51778row = 2399col = 3751print('投影坐标 -> 经纬度：')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]))