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
