GDAL笔记--chapter9
1.将独立的栅格波段合成为一张图像
#将独立的栅格波段合成一张图像 import os import numpy as np from osgeo import gdal os.chdir(r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Washington') band1_fn = r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Landsat\Washington\p047r027_7t20000730_z10_nn10.tif' band2_fn = r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Landsat\Washington\p047r027_7t20000730_z10_nn20.tif' band3_fn = r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Landsat\Washington\p047r027_7t20000730_z10_nn30.tif' in_ds = gdal.Open(band1_fn) in_band = in_ds.GetRasterBand(1) gtiff_driver = gdal.GetDriverByName('GTiff') out_ds = gtiff_driver.Create('nat_color.tif', in_band.XSize, in_band.YSize, 3, in_band.DataType) out_ds.SetProjection(in_ds.GetProjection()) print(in_ds.GetGeoTransform()) out_ds.SetGeoTransform(in_ds.GetGeoTransform()) in_data = in_band.ReadAsArray() out_band = out_ds.GetRasterBand(3) out_band.WriteArray(in_data) in_ds = gdal.Open(band2_fn) out_band = out_ds.GetRasterBand(2) out_band.WriteArray(in_ds.ReadAsArray()) #在数据集上调用ReadAsArray时,如果数据集有多个波段,则得到一个三维数组 #如果数据集是单波段,则返回二维数组。GetRasterBand的作用只是返回要查询的波段 in_ds = gdal.Open(band3_fn) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(in_ds.ReadAsArray()) out_ds.FlushCache() #确保数据写入磁盘,刷新缓存 for i in range(1, 4): out_ds.GetRasterBand(i).ComputeStatistics(False) #计算每个波段的统计数据 out_ds.BuildOverviews('average', [2,4,8,16,32]) #创建概视图 del out_ds #释放数据集
2.band.ReadAsArray(xoff,yoff,win_xsize,win_ysize,buf_xsize,buf_ysize,buf_obj)
#xoff是列起点、win_xsize是读取的列数、buf_xsize是输出数组的列数,如果和前者不同会重采样、buf_obj存放数组类型 #数据类型转换 data = np.empty((3, 6), dtype=float) band.ReadAsArrray(14000,6000,6,3, buf_obj=data) data = band.ReadAsArrray(1400,600,6,3).astype(float)
3.分块处理栅格
import os import numpy as np from osgeo import gdal os.chdir(r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Landsat\Washington') in_ds = gdal.Open('') in_band = in_ds.GetRasterBand(1) xsize = in_band.XSize ysize = in_band.YSize block_xsize, block_ysize = in_band.GetBlockSize() nodata = in_band.GetNoDataValue() #获取nodata out_ds = in_ds.GetDriver().Create('dem_feet.tif', xsize, ysize, 1, in_band.DataType) out_ds.SetProjection(in_ds.GetProjection()) out_ds.SetGeoTransform(in_ds.GetGeoTransform()) out_band = out_ds.GetRasterBand(1) for x in range(0, xsize, block_xsize): if x+block_xsize < xsize: cols = block_xsize else: cols = xsize-x for y in range(0, ysize, block_ysize): if y+block_ysize < ysize: rows = block_ysize else: rows = ysize-y data = in_band.ReadAsArray(x, y, cols, rows) data = np.where(data == nodata, nodata, data*3.28084) out_band.WriteArray(data, x, y) out_band.FlushCache() out_band.SetNoDataValue(nodata) #将nodata排除在外 out_band.ComputeStatistics(False) out_ds.BuildOverviews('average', [2, 4, 8, 16, 32]) del out_ds
4.GeoTransform
#使用现实世界的坐标 #geotransform #0、3 原点x\y坐标 现实坐标,一般是投影坐标 #1、5 像素的宽度和高度(高度负值) #2、4 x\y旋转 gt = ds.GetGeoTransform() #正变换,从图像坐标到现实坐标 inv_gt = gdal.InvGeoTransform(gt) #逆变换 offsets = gdal.ApplyGeoTransfrom(inv_gt, 465200, 5296000) #ApplyGeoTransform返回浮点数,若要传递给ReadAsArray需要转为整数 xoff, yoff = map(int, offsets) #python内置map函数,对数据做映射 value = band.ReadAsArrray(xoff, yoff, 1, 1)[0, 0] data = band.ReadAsArray() xoff, yoff = map(int, gdal.ApplyGeoTransform(inv_gt, 465200, 5296000)) value = data[yoff, xoff] #numpy数组[行,列],与gdal不同
5.保存图片的子集
#提取并保存图片的子集 import os from osgeo import gdal vashon_ulx, vashon_uly = 532000, 5262600 vashon_lrx, vashon_lry = 548500, 5241500 os.chdir(r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Landsat\Washington') in_ds = gdal.Open('nat_color.tif') in_gt = in_ds.GetGeoTransform() inv_gt = gdal.InvGeoTransform(in_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('Inverse geotransform failed') offsets_ul = gdal.ApplyGeoTransform(inv_gt, vashon_ulx, vashon_uly) offsets_lr = gdal.ApplyGeoTransform(inv_gt, vashon_lrx, vashon_lry) off_ulx, off_uly = map(int, offsets_ul) #取整得到行列值 off_lrx, off_lry = map(int, offsets_lr) rows = off_lry - off_uly columns = off_lrx - off_ulx gtiff_driver = gdal.GetDriverByName('GTiff') out_ds = gtiff_driver.Create('vashon2.tif', columns, rows, 3) out_ds.SetProjection(in_ds.GetProjection()) subset_ulx, subset_uly = gdal.ApplyGeoTransform(in_gt, off_ulx, off_uly) #这里作正变换是防止原左上角坐标落入像素中间的某个地方 print(subset_ulx) #计算出来的531995.25\5262624.75是原左上角坐标,正变换后不会落入像素中间。 print(subset_uly) out_gt = list(in_gt) #元组转列表,可修改 out_gt[0] = subset_ulx out_gt[3] = subset_uly out_ds.SetGeoTransform(out_gt) for i in range(1, 4): in_band = in_ds.GetRasterBand(i) out_band = out_ds.GetRasterBand(i) data = in_band.ReadAsArray(off_ulx, off_uly, columns, rows) out_band.WriteArray(data) out_ds.FlushCache() del out_ds
6.图像重采样
import os import numpy as np from osgeo import gdal os.chdir(r'') in_ds = gdal.Open(r'') in_band = in_ds.GetRasterBand(1) out_rows = in_band.YSize*2 out_columns = in_band.XSize*2 gtiff_driver = gdal.GetDriverByName('GTiff') out_ds = gtiff_driver.Create('band_resampled.tif', out_columns, out_rows) out_ds.SetProjection(in_ds.GetProjection()) geotransform = list(in_ds.GetGeoTransform()) geotransform[1] /= 2 geotransform[5] /= 2 out_ds.SetGeoTransform(geotransform) data = in_band.ReadAsArray(buf_xsize=out_columns, buf_ysize=out_rows) #指定buf_xsize缓冲,读取时会最邻近插值 out_band = out_ds.GetRasterBand(1) out_band.WriteArray(data) out_band.FlushCache() out_band.ComputeStatistics(False) out_ds.BuildOverviews('average', [2, 4, 8, 16, 32, 64]) del out_ds #重采样为更大的元素 data = np.empty((2, 3), np.int) #2行3列 band.ReadAsArray(1400, 6000, 6, 4, buf_obj=data)
7.ReadRaster\WriteRaster读取存储为字节序列
#ReadRaster\WriteRaster读取存储为字节序列 ds.ReadRaster(1400, 6000, 2, 2, band_list=[1]) ds.WriteRaster(1400, 6000, 6, 4, data, band_list=[1]) #字节序列ReadRaster import os import numpy as np from osgeo import gdal in_ds = gdal.Open(r'nat_color.tif') out_rows = int(in_ds.RasterYSize/2) #ds直接读取行列数 out_columns = int(in_ds.RasterXSize/2) num_bands = in_ds.RasterCount #ds读取波段数目 gtiff_driver = gdal.GetDriverByName('GTiff') out_ds = gtiff_driver.Create('nat_color_resampled.tif', out_columns, out_rows, num_bands) out_ds.SetProjection(in_ds.GetProjection()) geotransform = list(in_ds.GetGeoTransform()) geotransform[1] *= 2 geotransform[5] *= 2 out_ds.SetGeoTransform(geotransform) data = in_ds.ReadRaster(buf_xsize=out_columns, buf_ysize=out_rows) #利用较小的缓冲来读写数据,重采样 out_ds.WriteRaster(0, 0, out_columns, out_rows, data) #为输出数据源写入数据 out_ds.FlushCache() for i in range(num_bands): out_ds.GetRasterBand(i+1).ComputeStatistics(False) out_ds.BuildOverviews('average', [2, 4, 8, 16]) del out_ds
8.子数据集HDF(modis)
#子数据集HDF,GDAL的版本需要支持HDF from osgeo import gdal ds = gdal.Open(r'E:\桌面文件保存路径\gdal\osgeopy-data\osgeopy-data\Modis\MYD13Q1.A2014313.h20v11.005.2014330092746.hdf') subdatasets = ds.GetSubDatasets() #返回子数据集列表 print('Number of subdatasets:{}.format(len(subdatasets))') for sd in subdatasets: print('Name:{0}\nDescription:{1}\n'.format(*sd)) ndvi_ds = gdal.Open(subdatasets[0][0]) #第一个[0]读取第一个数据集,第二个则读取当前数据集的名称
9.xml保存影像(driver.CreateCopy)
ds = gdal.Open('listing9—6.xml') gdal.GetDriverByName('PNG').CreateCopy('liberty.png', ds) #driver.CreateCopy('', ds)用来复制ds