读取遥感影像数据并利用从参考影像的投影信息进行重投影,他们的空间范围一直
导入相关库
from osgeo import gdal
import numpy as np
import os
import fnmatch
import pandas as pd
import geopandas as gpd
import rasterio
from rasterstats import point_query
读写影像
def readData(envifile):
envidata = gdal.Open(envifile)
if envidata == None:
print(envifile+' Open Failed')
else:
print(envifile+' Open Success')
return envidata
def writeTiff(im_data, im_geotrans, im_proj, path):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float64
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
elif len(im_data.shape) == 2:
im_data = np.array([im_data])
im_bands, im_height, im_width = im_data.shape
# 创建文件
# 创建文件
driver = gdal.GetDriverByName("Gtiff")
dataset = driver.Create(path, int(im_width), int(
im_height), int(im_bands), datatype)
dataset.SetGeoTransform(im_geotrans) # 写入放射变换函数
dataset.SetProjection(im_proj) # 写入投影
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
print(path + 'Write Success!')
del dataset
重投影,根据一幅已有正确投影信息的影像对其他影像进行重投影
def reproject(inpath, outpath, referenceImg):
filelist = [] # 用来存所有文件
for root, dirs, files in os.walk(inpath):
for name in files:
if fnmatch.fnmatch(name, '*.tif'):
filelist.append(root + '\\' + name)
for imgfile in filelist:
# imgfile = 'H:\\CHL\\new\\01.tif'
dataset = readData(referenceImg)
Img_width = dataset.RasterXSize # 栅格矩阵的列数
Img_height = dataset.RasterYSize
Img_geotrans = dataset.GetGeoTransform() # 获取放射矩阵信息
Img_proj = dataset.GetProjection() # 获取投影信息
SavePath = '输出文件名.tif'
godataset = readData(imgfile)
data = godataset.ReadAsArray()
data[data<-99] = 0 #把小于-99的值设置为0
writeTiff(data, Img_geotrans, Img_proj, SavePath)
主函数
if __name__ == '__main__':
inpath = r'M:\01'
outpath = r'M:\11'
referenceImg = r'reference.tif'#有正确投影信息的影像
reproject(inpath, outpath, referenceImg)