GDAL笔记--chapter10

1.为栅格数据添加地面控制点

 

import shutil
from osgeo import gdal
orig_fn = r''
shutil.copy(orig_fn, fn)   #因为要更新，所以需要对文件做个备份
ds = gdal.Open(fn, gdal.GA_Update)
sr = osr.SpatialReference()
sr = osr.SetWellKnownGeogCS('WGS84')  #WGS84基准
gcps = [gdal.GCP(-111.93, 41.74, 0, 1078, 648),
        gdal.GCP(-111.90, 41.75, 0, 3531, 295)]

ds.SetGCPS(gcps, sr.ExportToWkt())  #将地面控制点附加到数据集
ds = None

#如果要使用一阶变换来创建地理变换，将地面控制点列表传递给GCPsToGeotransform,确保在数据集上设置了地理变换和投影信息

ds.SetProjection(sr.ExportToWkt())

ds.SetGeoTransform(gdal.GCPsToGeotransform(gcps)) #从地面控制点中创建地理变换，并设置在数据集上

 

2.栅格图像镶嵌

def get_extent(fn):
    ds = gdal.Open(fn)
    gt = ds.GetGeoTransform()
    return (gt[0], gt[3], gt[0]+gt[1]*ds.RasterXSize,
            gt[3]+gt[5]*ds.GetRasterYsize)

#将多幅影像镶嵌在一起
#Transform类可以帮助在真实世界坐标和图像坐标之间转换
#还可以在两个不同的栅格之间进行像素偏移
import os
import glob
import math
os.chdir(r'')
in_files = glob.glob('o*.tif')
min_x, max_y, max_x, min_y = get_extent(in_files[0])
for fn in in_files[1:]:
    minx, maxy, maxx, miny = get_extent(fn)
    min_x = min(min_x, minx)
    max_y = max(max_y, maxy)
    max_x = max(max_x, maxx)
    min_y = min(min_y, miny)

in_ds = gdal.Open(in_files[0])
gt = in_ds.GetGeoTransform()  #这里每幅图的分辨率都是一样的
rows = math.ceil((max_y-min_y)/-gt[5]) #计算镶嵌后的行列数
columns = math.ceil((max_x-min_x)/gt[1])

driver = gdal.GetDriverByName('gtiff')
out_ds = driver.Create('mosaic.tif', columns, rows)
out_ds.SetProjection(in_ds.GetProjection())
out_band = out_ds.GetRasterBand(1)

gt = list(in_ds.GetGeoTransform())
gt[0], gt[3] = min_x, max_y  #设置新的起始点坐标
out_ds.SetGeoTransform(gt)

for fn in  in_files:
    in_ds = gdal.Open(fn)
    trans = gdal.Transformer(in_ds, out_ds, []) #第三个参数设置转换器选项，这里不需要
    #False表示计算从源栅格到目标栅格的偏移
    success, xyz = trans.TransformPoint(False, 0, 0)  #这里转换的是图像的左上角行列号坐标，大图将会从x,y处开始读取数据，所以输入应该是(0,0)
    x, y, z = map(int, xyz)
    data = in_ds.GetRasterBand(1).ReadAsArray()
    out_band.WriteArray(data, x, y)

del in_ds, out_band, out_ds

#将图像坐标转换成现实世界坐标，这里out_ds是输入吧？
trans = gdal.Transformer(out_ds, None, [])
success, xyz = trans.TransformPoint(0, 1078, 648)

3.为栅格数据添加颜色映射

#为栅格数据添加颜色映射
import os
from osgeo import gdal
os.chdir('E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Switzerland')
original_ds = gdal.Open('dem_class.tif')
driver = gdal.GetDriverByName('gtiff')
ds = driver.CreateCopy('dem_class2.tif', original_ds)
#备份文件
band = ds.GetRasterBand(1)
colors = gdal.ColorTable()
colors.SetColorEntry(1, (112, 153, 89))  #第一个参数是value，第二个是RGB分量
colors.SetColorEntry(2, (242, 238, 162))
colors.SetColorEntry(3, (242, 206, 133))
colors.SetColorEntry(4, (194, 140, 124))
colors.SetColorEntry(5, (214, 193, 156))

band.SetRasterColorTable(colors) 
#把颜色映射添加到波段
band.SetRasterColorInterpretation(gdal.FCI_PaletteIndex)
#将颜色解译为调色板
del band, ds

#颜色表的更改
ds = gdal.Open('')
band = ds.GetRasterBand(1)
colors = band.GetRasterColorTable()  #获取颜色表
colors.SetColorEntry(5, (250, 250, 250))  #更改颜色
band.SetRasterColorTable(colors)  #添加到颜色表
del band, ds


#alpha代表不透明度，在创建时需要指定第二个波段是一个alpha波段
#第一个波段存储像素值，第二个是透明度
ds = driver.Create('dem_class4.tif', original_ds.RasterXSize, 
                    original_ds.RasterYSize, 2 , gdal.GDT_Byte, ['ALPHA=YES'])
#像素为5透明度65，其余为255，写入band2，最后解译透明度波段
import numpy as np
data = band1.ReadAsArray()
data = np.where(data==5, 65, 255)
band2.WriteAsArray(data)
band2.SetRasterColorInterpretation(gdal.GCI_AlphaBand)

4.GetHistogram函数

 

import os
from osgeo import gdal
os.chdir('E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Switzerland')
ds = gdal.Open('dem_class2.tif')
band = ds.GetRasterBand(1)
approximate_hist = band.GetHistogram()
exact_hist = band.GetHistogram(approx_ok=False)
print('Approximate:', approximate_hist[:7], sum(approximate_hist))
print('Exact:', exact_hist[:7], sum(exact_hist))

hist = band.GetHistogtam(0.5, 6.5, 3, approx_ok=False)
band.SetDefaultHistogram(1, 6, hist)  #用最小像素值、最大像素值和计数列表设置默认直方图
min_val, max_val, n, hist = band.GetDefaultHistogram() #返回最小像素、最大像素、组距数、一个计数列表的元组
print(hist)

 

5.为栅格数据添加属性表

 

#为栅格数据添加属性表
os.chdir('')
ds = gdal.Open('dem_class2.tif')
band = ds.GetRasterBand(1)
band.SetNoDataValue(-1)  #不存在的值设置为nodata

rat = gdal.RasterAttributeTable()   #创建属性表
rat.CreateColumn(                  #创建列
    'Value', gdal.GFT_Integer, gdal.GFU_Name     #表示含像素值的列
)
rat.CreateColumn(
    'Count', gdal.GFT_Integer, gdal.GFU_PixelCount  #表示直方图计数
)
rat.CreateColumn(
    'Elevation', gdal.GFT_String, gdal.GFU_Generic  #表示描述信息
)
rat.SetRowCount(6)   #创建行

rat.WriteArray(range(6), 0)   #第一列，012345
rat.WriteArray(band.GetHistogram(-0.5, 5.5, 6, False, False), 1)  #第二列，直方图计数Count
rat.SetValueAsString(1, 2, '0-800')   #第三列elevation
rat.SetValueAsString(2, 2, '800-1300')
rat.SetValueAsString(3, 2, '1300-2000')
rat.SetValueAsString(4, 2, '2000-2600')
rat.SetValueAsString(5, 2, '2600+')

band.SetDefaultRAT(rat)  #添加属性表到波段
band.SetNoDataValue(0)  #恢复nodata
del band, ds

 

6.使用xml定义有一个波段的VRT数据集

 

<VRTDataset rasterXSize='8849' rasterYSize='8023'>
    <SRS>
        PROJCS['WGS 84 /UTM zone 10N', GEOGCS['WGS 84',DATUM['WGS_1984'...]
    </SRS>
    <GeoTransform>343724.25,28.5,0,5369585.25,0,-28.5</Geotransform>
    <VRTRasterBand dataType='Byte' band='1'>
        <SimpleSource>
            <SourceFileName relativeToVRT='1'>
                nat_color.tif
            </SourceFilename>
            <SourceBand3>3</SourceBand>
            <SourceProperties RasterXSize='8849' RasterYSize='8023'
                DataType='Byte' BlockXSize='8849' BlockYSize='1'/>
        </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

 

7.VRT创建

#VRT
from osgeo import gdal
import os
os.chdir('')
tmp_ds = gdal.Open('')
driver = gdal.GetDriverByName('vrt')
ds = driver.Create('nat_color.vrt', tmp_ds.RasterXSize, tmp_ds.RasterYSize, 3)
ds.SetProjection(tmp_ds.GetProjection())
ds.SetGeoTransform(tmp_ds.GetGeoTransform())

#将链接添加到栅格。创建键值source_0,包含文件名的xml字符串
metadata = {'source_0':xml.format('....tif')}
ds.GetRasterBand(1).SetMetadata(metadata, 'vrt_sources')

metadata = {'source_0':xml.format('....tif')}
ds.GetRasterBand(2).SetMetadata(metadata, 'vrt_sources')

metadata = {'source_0':xml.format('....tif')}
ds.GetRasterBand(3).SetMetadata(metadata, 'vrt_sources')

8.VRT裁剪影像

#VRT影像裁剪
import os
from osgeo import gdal

os.chdir('')
tmp_ds = gdal.Open('')
tmp_gt = tmp_ds.GetGeoTransform()
inv_gt = gdal.InvGeoTransform(tmp_gt)
if gdal.VersionInfo()[0]=='1':
    if inv_gt[0]==1:
        inv_gt = inv_gt[1]
    else:
        raise RuntimeError('Inverse geotransform failed')
elif inv_gt is None:
    raise RuntimeError('failed')

vashon_ul = (532000, 5262600)
vashon_lr = (548500, 5241500)
ulx, uly = map(int, gdal.ApplyGeoTransform(inv_gt, *vashon_ul))   #计算裁剪部分左上角图像坐标
lrx, lry = map(int, gdal.ApplyGeoTransform(inv_gt, *vashon_lr))   #裁剪部分右下角
rows = lry - uly
columns = lrx - ulx
gt = list(tmp_gt)
gt[0] += gt[1]*ulx  #原gt[0]是原图的左上角，计算后为裁剪的左上角
gt[3] += gt[5]*uly  #裁剪的右下角

ds = gdal.GetDriverByName('vrt').Create('vashon.vrt', columns, rows, 3)
ds.SetProjection(tmp_ds.GetProjection())
ds.SetGeoTransform(tmp_ds.GetGeoTransform())  #地理变换为裁剪对象的地理变换

xml=''             #描述波段的xml文件
<SimpleSource>
    <SourceFilename relativeToVRT='1'>{fn}</SourceFilename>
    <SourceBand>{band}</SourceBand>
    <SrcRect xOff='{xoff}' yOff='{yoff}'
            xSize='{cols}' ySize='{rows}'/>
    <DstRect xOff='0' yOff='0'
            xSize='{cols}' ySize='{rows}'/>
</SimpleSource>

data = {'fn':'nat_color.tif', 'band':1,   #将数据插入到xml文件
        'xoff':ulx, 'yoff':uly,
        'cols':columns, 'rows':rows}

meta = {'source_0': xml.format(**data)}   #波段1插入到VRT
ds.GetRasterBand(1).SetMetadata(meta, 'vrt_sources')

data['band'] = 2
meta = {'source_0': xml.format(**data)}
ds.GetRasterBand(2).SetMetadata(meta, 'vrt_sources')

data['band'] = 3
meta = {'source_0': xml.format(**data)}
ds.GetRasterBand(3).data(meta, 'vrt_sources')

del ds, tmp_ds

9.创建问题格式

#创建问题格式
ds = gdal.Open('vashon.vrt')
gdal.GetDriverByName('jpeg').CreateCopy('vashon.jpg', ds)

10.栅格数据重投影

#栅格数据的像元易弯曲和移动，通常使用最近领域插值法、双线性插值和三次卷积插值
#方法1.gdalwarp 2.AutoCreateWarpedVRT
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
#目标srs是WGS84基准的地理坐标系，即未投影。将把UTM投影转换为地理坐标系
old_ds = gdal.Open('nat_color.tif')
vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(),
    gdal.GRA_Bilinear)
#没有在磁盘上创建VRT文件，而是返回一个数据集对象，可以用CreateCopy保存为其他格式
gdal.GetDriverByName('gtiff').CreateCopy('nat_color_wgs84.tif', vrt_ds)

11.回调函数

import sys
def my_progress(complete, message, progressArg=0.02):
    if not hasattr(my_progress, 'last_progress'):
        sys.stdout.write(message)
        my_progress.last_progress=0
    if complete >= 1:
        sys.stdout.write('done/n')
        del my_progress.last_progress
    else:
        progress = divmod(complete, progressArg)[0]
        while my_progress.last_progress < progress:
            sys.stdout.write('.')
            sys.stdout.flush()
            my_progress.last_progress += 1