python GIS 影像处理

打开栅格并获取信息

from osgeo import gdal
fn = “E:\landuse_biodiversity\Data\Biome.tif”

ds = gdal.Open(fn)
print(ds.RasterXSize, ds.RasterYsize)
print(ds.GetProjection())
print(ds.GetGeoTransform()) #地理变换，【左上角的精度，X分辨率，X旋转角度，左上角的纬度，Y旋转角度，Y分辨率】
print(ds.RasterCount())

band1 = ds.GetRasterBand(1)
print(band1.NoDataValue())
print(band1.GetMinimum())
print(band1.GetMaximum())
print(band1.GetUnitType())

读取栅格数据为数组

from osgeo import gdal

fn="c:/temp/test.tif"
ds = gdal.Open(fn)

band1= ds. GetRasterBand(1).ReadAsArray()
print(band1.shape)
print(band1)

创建栅格数据

from osgeo import gdal

fn = "c:/temp/new_raster.tif"
fn2 = "C:/temp/raster.tif"

ds = gdal.Open(fn2)
geot = ds.GetGeoProjection()
proj = ds.GetProjection()

dirver_tiff  = gdal.GetDriverByName("GTiff")
new_ds=driver_tiff.Create(fn , xsize=ds.RasterXSize, ysize=ds.RasterYSize, bands=1, etype=gdal.GDT_Float32)
new_ds.SetGeoTransform(geot)
new_ds.SetProjection(proj)

复制一个栅格

from osgeo import gdal

fn = "dem.tif"
fn_new = "dem_copy" 

ds = gdal.Open(fn)

driver_tiff = gdal.GetDriverByName("GTiff")

ds_copy = driver_tiff.CreateCopy(fn_new,ds , strict=0)
ds_copy = None
#iface.addRasterLayer(fn_new)#将fn_new在qgis中打开

改变栅格的值

from osgeo import gdal
fn = "C:/temp/test_raster.tif"
fn_template = "C:/temp/elevation/ned10m.tif"

ds_temp = gdal.Open(fn_template)

driver_tiff = gdal.GetDriverByName("GTiff")

ds = driver_tiff.Create(fn, xsize=10, ysize=10, bands=1, eType=gdal.GDT_Float32)

ds.SetGeoTransform(ds_temp.GetGeoTransform())
ds.SetProjection(ds.temp.GetProjection())

band1 = ds.GetRasterband(1).ReadAsArray()

band1[0][0] = 100.0
band1[5][5] = 50.0
band1[1][0] = 10.0

ds.GetRasterBand(1).WriteArray(band1)

#iface.addRasterLayer(fn_new)#将fn_new在qgis中打开

栅格数据的No data

from osgeo import gdal
import numpy as np

fn = "dem.tif"
fn_new = "dem_copy" 

driver = gdal.GetDriverByName("GTfiff")
copyds = gdal.Open(fn)
ds = driver.CreateCopy(fn_new, copyds, strict=0)

ds = gdal.Open(fn_new)
ndv = ds.GetRasterBand(1).GetNoDataValue()
band1 = ds.GetRasterBand(1).ReadAsArray()
band1 = band1*2.0
print(band1)
band1 = ds.GetRasterBand(1).ReadAsArray()
band1[band1==ndv] = np.nan
band1 = band1*2.0
ds.GetRasterband(1).WriteArray(band1)
print(band1)