Python GIS 数据处理

最新推荐文章于 2024-09-15 22:31:42 发布

NeverDelay001

最新推荐文章于 2024-09-15 22:31:42 发布

阅读量354

点赞数 6

文章标签： python

本文链接：https://blog.csdn.net/weixin_42622230/article/details/135713802

版权

本文介绍了如何使用GDAL库在Python中进行栅格图层相减、计算NDVI（归一化差值指数）以及进行地理空间数据的投影转换和裁剪。通过实例展示了如何打开和处理遥感影像数据，如DEM和Landsat数据。

摘要由CSDN通过智能技术生成

栅格图层相减

from osgeo import gdal


def openRaster(fn, access=0):
	ds = gdal.Open(fn,access)
	if ds is None:
		print("Error opening raster dataset")
	return ds
def getRasterBand(ds, band=1, access=0):
	ds = openRaster(fn, access)
	band = ds.GetRasterBand(1).ReadAsArray()
	return band

def creatRasterFromCopy(fn, ds, data, driverFmt="GTiff"):
	dirver = gdal.GetDriverByName(diverFmt)
	outds = driver.CreateCopy(fn, ds, strict=0)
	outds.GetRasterBand(1).WriteArray(data)
	ds = None
	outds = None

def subtractRasters(fn1 , fn2)
	dat1 = getRasterBand(fn1)
	dat2 = getRasterBand(fn2)
	datout = dat2 - dat1
	createRasterFromCopy(fnout, gdal.open(fn1), datout)

fn1 = ”dem1.tif“
fn2 = "dem2.tif"
fnout = "subtract.tif"
subtractRasters(fn1, fn2 fnout)
#ifcace.addRasterLayer(fnout)

计算NDVI

from osgeo import gdal
import numpy as np

def openRaster(fn, access=0):
	ds = gdal.Open(fn,access)
	if ds is None:
		print("Error opening raster dataset")
	return ds
	
def getRasterBand(ds, band=1, access=0):
	ds = openRaster(fn, access)
	band = ds.GetRasterBand(1).ReadAsArray()
	return band

def creatRasterFromCopy(fn, ds, data, driverFmt="GTiff"):
	dirver = gdal.GetDriverByName(diverFmt)
	outds = driver.CreateCopy(fn, ds, strict=0)
	outds.GetRasterBand(1).WriteArray(data)
	ds = None
	outds = None

def creatRasterFromTemplate(fn, ds, data, ndv=-9999.0, driverFmt="GTiff"):
	driver = gdal.GetDriverByName(diverFmt)
	outds = driver.Create(fn, xsize=ds.RasterXSize, ysize=ds.RasterYSize, band=1, eType=gdal.GDT_Float32)
	outds.SetGeoTransform(ds.GetGeoTransform())
	outds.SetProjection(ds.GetProjection())
	outds.GetRasterBand(1).SetNoDataValue(ndv)
	outds.GetRasterBand(1).WriteArray(data)
	outds = None
	ds = None

def ndvi(nirband, redband, ndv=-9999.0):
	nirband[nirband < 0] = np.nan
	nirband[nirband > 10000] = np.nan
	redband[redband < 0] = np.nan
	redband[redband > 10000] = np.nan
	ndviband = (nirband-redband)/(nirband+redband)
	ndviband[np.isnan(ndivband)=ndv
	return ndviband 

infn = "landsat7_test.tif"
outfn = "ndvi.tif"
redband = getRasterBand(infn, 4)
nirband = getRasterBand(infn, 5)
ndviband = ndvi(nirband, redband)
createRasterFromTemplate(outfn, gdal.Open(infn), ndviband)

投影转换

from osgeo import gdal

fn = "ned10m.tif"
out = "ned10m_4326.tif"

gdal.Warp(out, fn, dstSRS="EPSG:4326",xRes=0.1, yRes=0.1)#width=1065,height=426

裁剪栅格影像

from osgeo import gdal

fn = "ned10m.tif"
fn_clip = "ned10m_clip"

# 左上角坐标不变化
ds = gdal.Open(fn)
newRows = 700  #y coordinate
newCols = 500  #x coordinate

geot = ds.GetGeoTransform()

driver = gdal.GetDriverByName("GTiff")
ds_clip = driver.Create(fn_clip, xsize=newCols, ysize=newRows,bands=1,eType=gdal.GDT_Float32)

ds_clip.SetGeoTransform(geot)

#data = ds.GetRasterBand(1).ReadAsArray()
#ds_clip.GetRasterBand(1).WriteArray(data[:newRows, :newCols])
#ds_clip.GetRasterBand(1).SetNoDataValue(ds.GetRasterBand(1).GetNoDataValue())

#例子2，裁剪时左上角坐标变化
xOffset = 2500
yOffset = 3000

geoNew = geot

geotNew[0] = geotNew[0] + xOffset
geotNew[3] = geotNew[3] + yOffset

ds_clip.SetGeoTransform(geotNew)
startRow = round(yOffset/abs(geot[5]))
startCol = round(xOffset/abs(geot[1]))

data = ds.GetRasterBand(1).ReadAsArray()
ds_clip.GetRasterBand(1).WriteArray(data[startRow:startRow+newRows,startCol:startCol+newCols])
#iface.addRasterLayer(fn_clip)