Python GDAL学习笔记(二)

一、自定义空间坐标系

使用EPSG或EPSGA编号进行初始化

EPSG：欧洲石油勘探组织（European Petroleum Survey Group，EPSG)，EPSG编译并传播了EPSG大地参数集、广泛使用的地球椭球体数据库、大地基准、地理和投影坐标系、测量单位等。该数据集是坐标参考系统和坐标转换定义的集合，这些定义可以应用在全球、区域、国家或局部地区。目前，这个数据集被广泛接受并使用。EPSG通过一个Web发布平台http://www.epsg.org/对该数据集进行分发。在该数据库中存储的椭球体，投影坐标系等不同组合都对应着不同的ID号，这个号在EPSG中被称为EPSGcode，它代表特定的椭球体、单位、地理坐标系或投影坐标系等信息。

from osgeo import gdal, osr

srs = osr.SpatialReference()
srs.ImportFromEPSG(4509) #WGS 84，其编码为4326
print(srs)

PROJCS["CGCS2000 / Gauss-Kruger CM 117E",
    GEOGCS["China Geodetic Coordinate System 2000",
        DATUM["China_2000",
            SPHEROID["CGCS2000",6378137,298.257222101,
                AUTHORITY["EPSG","1024"]],
            AUTHORITY["EPSG","1043"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4490"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Northing",NORTH],
    AXIS["Easting",EAST],
    AUTHORITY["EPSG","4509"]]

使用WKT字符串进行初始化

首先获取一幅已经进行投影的proj的Wkt信息。

from osgeo import gdal, osr

filename = r'F:\GDAL learning\Landsat8.tif'
dataset = gdal.Open(filename, gdal.GA_ReadOnly)
proj = dataset.GetProjection()
print(proj)

PROJCS["unnamed",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","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

进行自定义参考系统

srs = osr.SpatialReference()
srStr = 'PROJCS["unnamed",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","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
srs.ImportFromWkt(srStr)
print(srs)

PROJCS["unnamed",
    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","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["latitude_of_center",0],
    PARAMETER["longitude_of_center",105],
    PARAMETER["standard_parallel_1",25],
    PARAMETER["standard_parallel_2",47],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

使用PROJ字符串进行初始化

from osgeo import gdal, osr

filename = r'F:\GDAL learning\Landsat8.tif'
dataset = gdal.Open(filename, gdal.GA_ReadOnly)
proj = dataset.GetProjection()

srs = osr.SpatialReference()
srs.ImportFromWkt(proj)
print(srs.ExportToproj4())

+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

srs = osr.SpatialReference()
srs.ImportFromProj4('+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 
                    +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
print(srs)

参数	描述
+a	椭球体的长半轴
+axis	坐标轴方向
+b	椭球体的短半轴
+ellps	椭球体的名称
+k	缩放比例（新版本中该参数已经被弃用）
+k_0	缩放比例
+lat_0	起始纬度
+lon_0	中央经线
+lon_wrap	区域中央经线的范围
+over	是否允许地图的范围不在-180度到180度之间。
+pm	指定本初子午线，用以说明所使用的坐标系中的本初子午线与格林威治本初子午线之间的偏移，一般是一个城市的名字。我们使用greenwich就行了。
+proj	投影名称
+units	单位，米、英尺等
+vunits	垂直方向的单位.
+x_0	伪东向
+y_0	伪北向

附录

osr.SpatialReference()定义一个Spref对象
srs.ImportFromEPSG(4326)通过EPSG编码为spref添加坐标信息

proj = dataset.GetProjection() 获取栅格Wkt投影信息
srs.ImportFromWkt(proj)通过Wkt为Spref添加投影信息

srs.ExportToProj4()空间参考输出为proj4格式
srs.ImportFromProj4('+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')通过Proj4为Spref添加投影信息

二、遥感影像几何校正

利用卫星自带的地理定位文件进行几何校正

查看数据子数据集

from osgeo import gdal, osr

Datafilename = 'F:\MOD11__\MOD11_L2.A2015061.0250.006.2016215171955.HDF'

Geodataset = gdal.Open(Geofilename)
dataset = gdal.Open(Datafilename)

subdatasets = dataset.GetSubDatasets()
SubDatasetsNum =  len(subdatasets)

print("子数据集一共有{0}个: ".format(SubDatasetsNum))
for i in range(SubDatasetsNum):
    print(subdatasets[i])

子数据集一共有16个: 
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:LST', '[2030x1354] LST MOD_Swath_LST (16-bit unsigned integer)')
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:QC', '[2030x1354] QC MOD_Swath_LST (16-bit unsigned integer)')
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:Error_LST', '[2030x1354] Error_LST MOD_Swath_LST (8-bit unsigned integer)')
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:Emis_31', '[2030x1354] Emis_31 MOD_Swath_LST (8-bit unsigned integer)')
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:Emis_32', '[2030x1354] Emis_32 MOD_Swath_LST (8-bit unsigned integer)')
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:View_angle', '[2030x1354] View_angle MOD_Swath_LST (8-bit unsigned integer)')
('HDF4_EOS:EOS_SWATH:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":MOD_Swath_LST:View_time', '[2030x1354] View_time MOD_Swath_LST (8-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":0', '[406x271] Latitude (32-bit floating-point)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":1', '[406x271] Longitude (32-bit floating-point)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":2', '[2030x1354] LST (16-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":3', '[2030x1354] QC (16-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":4', '[2030x1354] Error_LST (8-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":5', '[2030x1354] Emis_31 (8-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":6', '[2030x1354] Emis_32 (8-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":7', '[2030x1354] View_angle (8-bit unsigned integer)')
('HDF4_SDS:UNKNOWN:"F:\\MOD11__\\MOD11_L2.A2015061.0250.006.2016215171955.HDF":8', '[2030x1354] View_time (8-bit unsigned integer)')

进行几何校正

subName = subdatasets[0][0] # 代表第一个子数据集数据 LST
subds = gdal.Open(subName)
#array = subds.ReadAsArray()
correctdata = gdal.Warp("F:\GDAL learning\mod11_geo_correct.tif",  
          				subds, dstSRS="EPSG:4326" ,format='GTiff', geoloc=True)
del correctdata

校正结果

速度可以，但是没有去除bowtie效应。

使用GCPs校正

from osgeo import gdal, osr

filename = 'F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_0250M_MS.HDF'

dataset = gdal.Open(filename)
SubDatasets = dataset.GetSubDatasets()
SubDatasetNum = len(SubDatasets)

print('子数据集一共{i}个：'.format(i = SubDatasetNum))
for i in range(SubDatasetNum):
    print(SubDatasets[i])
    
Metadata = dataset.GetMetadata()
MetadataNum = len(Metadata)
print('元数据一共有{0}个: '.format(MetadataNum))
for key,value in Metadata.items():
    print('{key}:{value}'.format(key = key, value = value))

Latitudes = Metadata['Orbit_Point_Latitude']
LatitudeList = Latitudes.split(' ')
Longitudes = Metadata['Orbit_Point_Longitude']
LongitudeList = Longitudes.split(' ')

SubdataStr = SubDatasets[9][0]
SubData = gdal.Open(SubdataStr)

row = SubData.RasterYSize
col = SubData.RasterXSize

#虽然前面获取了行列号和经纬度列表，但不确定对应关系，这里还是手动输了一下控制点。
#四个点有点少，应该从经纬度文件中每隔50个取一个控制点来校正会比较好。
#因为控制是经纬度对应的行列号，好像不可以转换投影坐标（不确定）。
Gcps = [gdal.GCP(140.75398, 38.429569, 0, 8191, 7999),
        gdal.GCP(109.42703, 33.700436, 0, 0, 7999),
        gdal.GCP(140.17953, 55.829273, 0, 8191, 0),
        gdal.GCP(98.286346, 49.495743, 0, 0, 0)]

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
SubData.SetGCPs(Gcps, srs.ExportToWkt())

dst_ds = gdal.Warp('MERSI_Band_Gcpgeo.tif', SubData, format='GTiff', tps = True,
                      xRes = 0.0025, yRes = 0.0025, srcNodata = 0, dstNodata = 0,
                resampleAlg = gdal.GRA_Bilinear, outputType=gdal.GDT_Float32)
del dst_ds                

校正结果

速度很快，一景只用了三分钟，ENVI建立查找表校正要花费三四个小时呢。
但因为只是单纯的几何校正，而且控制点取得很少，所以影像边缘bowtie很严。

三、空间参考系统转换

转换为Albers投影坐标

from osgeo import gdal, osr

filename = 'F:\GDAL learning\MERSI_Band_Gcpgeo.tif'
dataset = gdal.Open(filename)

srs = osr.SpatialReference()
srs.ImportFromProj4('+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')

repro_data = gdal.Warp('Albers.tif', dataset, dstSRS = srs, 
						xRes = 250,yRes = 250, 
						srcNodata = 0, dstNodata = 0,
                        resampleAlg = gdal.GRA_Bilinear,
                        outputType = gdal.GDT_Float32)
del repro_data

重投影结果

四、Warp参数详解

对于栅格数据的重投影可以使用Warp来实现
Warp(destNameOrDestDS, srcDSOrSrcDSTab,**kwargs)

• destNameOrDestDS是数据集变量或者新保存的目标文件地址。
• srcDSOrSrcDSTab可以是源数据集或文件名，即被重投影的数据集，也可以是包含了若干数据集或文件名的一个array。
• kwargs关键字参数是可以是一个WarpOptions对象，也可以是一个字符串数组，包含了WarpOptions的若干条件，具体内容与WarpOptions的参数一致。如果已经使用了一个WarpOptions对象，则所使用的其他关键字参数都会被忽略掉。
  W.arpOptions(options=[],
  format=None,
  outputBounds=None,
  outputBoundsSRS=None,
  xRes=None, yRes=None,
  targetAlignedPixels=False,
  width=0, height=0,
  srcSRS=None, dstSRS=None,
  srcAlpha=False, dstAlpha=False,
  warpOptions=None,
  errorThreshold=None,
  warpMemoryLimit=None,
  creationOptions=None,
  outputType=GDT_Unknown,
  workingType=GDT_Unknown,
  resampleAlg=None,
  srcNodata=None, dstNodata=None,
  multithread=False,
  tps=False,
  rpc=False,
  geoloc=False,
  polynomialOrder=None,
  transformerOptions=None,
  cutlineDSName=None,
  cutlineLayer=None,
  cutlineWhere=None,
  cutlineSQL=None,
  cutlineBlend=None,
  cropToCutline=False,
  copyMetadata=True,
  metadataConflictValue=None,
  setColorInterpretation=False,
  callback=None, callback_data=None
  )
  下面是对其参数的解释。
  options — 可以是一个字符串数组，一个字符串或者令其为空值，但是使用后面其他的参数来定义。
  format — 输出的格式 (例如"GTiff"等)。
  outputBounds — 在目标空间参考系统的输出数据集的范围，形式为 (minX,minY, maxX, maxY) 。
  outputBoundsSRS — 如果在dstSRS中没有定义的话，使用这个关键字定义输出数据集的边界的空间参考系统。
  xRes, yRes — 在目标参考系统中的像元大小。
  targetAlignedPixels —是否强制输出边界为输出分辨率的倍数。
  width — 输出栅格的像素列数。
  height — 输出栅格的像素行数。
  srcSRS —源空间参考系统。
  dstSRS — 输出空间参考系统。
  srcAlpha — 是否强制将输入数据集的最后一个波段作为alpha波段。
  dstAlpha — 是否强制创建一个输出数据集的alpha波段。
  outputType — 输出类型 (例如gdal.GDT_Byte等)
  workingType — working type (gdal.GDT_Byte, etc…)
  warpOptions —变形选项列表。
  errorThreshold --近似转换的误差阈值(用像素表示) 。
  warpMemoryLimit — 工作缓存大小，单位是bytes。
  resampleAlg — 重采样模式。
  creationOptions — 创建选项列表。
  srcNodata — 源数据的nodata值。
  dstNodata — 输出数据的nodata值。
  multithread — 是否多线程计算和输入输出操作。
  tps— 是否使用Thin Plate Spline GCP 转换器。
  rpc— 是否使用RPC转换器。
  geoloc — 是否使用GeoLocation数组转换器。
  polynomialOrder — 多项式GCP插值的阶数。
  transformerOptions — 转换参数
  cutlineDSName — 剪切线数据集名称。这里的剪切线是指对影像进行剪切的时候所使用的矢量图层。
  cutlineLayer — 剪切线图层名称。
  cutlineWhere — 剪切线的WHERE语句。
  cutlineSQL — 剪切线的SQL 语句。
  cutlineBlend — 以像素为单位的剪切线混合距离。
  cropToCutline — 是否使用剪切线的extent作为输出的界线。
  copyMetadata — 是否拷贝源数据的元数据。
  metadataConflictValue — 元数据冲突值。
  setColorInterpretation — 是否强制将输入波段的颜色解释赋予输出波段。
  callback — 回调函数。
  callback_data — 回调函数数据

resampling_method

GRA_NearestNeighbour = 0
GRA_Bilinear = 1
GRA_Cubic = 2
GRA_CubicSpline = 3
GRA_Lanczos = 4
GRA_Average = 5
GRA_Mode = 6
GRA_Gauss = 7
GRA_Max = 8
GRA_Min = 9
GRA_Med = 10
GRA_Q1 = 11
GRA_Q3 = 12