该博客详细介绍了如何使用GDAL库读取、查看和处理遥感数据，包括读取不同来源的遥感影像，如无人机、高分2号和北京2号数据，以及普通3波段图像。内容涵盖获取数据信息（如高度、宽度、波段数等）、将数据转换为numpy数组、分波段操作、数据类型转换以及数据保存为TIFF文件。此外，还展示了如何分块读取和写入数据，以及通过Counter统计数组元素频率。
gdal读取了遥感数据后就变成了numpy数组了，之后就是numpy的操作了，numpy操作了之后再加坐标保存就行了，实际上很简单，读取和写入看下面的就行了
1 图像操作
1.1 查看图像的信息
from osgeo import gdal
import gdal
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from collections import Counter
# 一、 读取数据
# 1 读取无人机
# 4波段
#r"J:\4影像数据大全\1无人机数据\巴南区无人机和高分影像\1巴南区无人机影像\正射影像\77-04.tif 就能读取
#r"J:\4影像数据大全\1无人机数据\巴南区无人机和高分影像\1巴南区无人机影像\正射影像\75-01.tif gdal2.3.3版本不能读取
# 3波段
# H:\01HTutorWork\3GF2\3Code\case\anli12_fx\fx163_1.tif
#in_ds = gdal.Open(r"J:\4影像数据大全\1无人机数据\巴南区无人机和高分影像\1巴南区无人机影像\正射影像\83-02.tif")
# 2 读取高分2号
#in_ds = gdal.Open(r"J:\4影像数据大全\3高分2号数据\2辽宁高分数据\2辽宁处理好的20景影像20210527\GF2_PMS1_E124.0_N41.9_20180924_L1A0003483014-MSS1.tif")
# 3 读取北京2号
#in_ds = gdal.Open(r"J:\4影像数据大全\2北京2号数据\1北京2号数据46景\8个县数据成果46景\TRIPLESAT_1_PSH_L4_20200415023916_002A77VI_016.pix")
# 4 读取普通的3波段jpg
#in_ds = gdal.Open(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\2arcgisl罗田县\2labelimg\lt497908_PolygonToRaster.tif")
#in_ds = gdal.Open(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\2arcgislLuotian\3imageDataset\1lt497908\labels\128.tif")
#in_ds = gdal.Open(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\1fujian_shuncang\5dataset16\2labelsPositive_N_png\p0lt497908n958813.png")
in_ds = gdal.Open(r"H:\fenge\shishi\ce255mian.tif")
# 二、 参看数据信息，获取numpy值
print("open tif file succeed")
height = in_ds.RasterYSize # 获取数据高度
print("高_行：",height)
width = in_ds.RasterXSize # 获取数据宽度
print("宽_列：",width)
outbandsize = in_ds.RasterCount # 获取数据波段数
print("波段数：",outbandsize)
im_geotrans = in_ds.GetGeoTransform() # 获取仿射矩阵信息
print("仿射信息：",im_geotrans)
im_proj = in_ds.GetProjection() # 获取投影信息
print("投影信息：",im_proj)
datatype = in_ds.GetRasterBand(1).DataType # 这个数据类型没搞清楚，1，2，3代表什么https://blog.youkuaiyun.com/u010579736/article/details/84594742
print("数据类型：",datatype)
# 2 分块或者分波段得到numpy数组
# img_array = in_ds.ReadAsArray(0, 0, width, height) # 将数据写成数组，分块读取，指定区域
img_array = in_ds.ReadAsArray() # 将数据写成数组，读取全部数据，结果是numpy数组
print("数组类型名",img_array.dtype.name)
print("数组形状",img_array.shape)
print("全部波段最大值：", np.max(img_array))
print("全部波段最小值", np.min(img_array))
print("全部波段唯一值", np.unique(img_array))
#print('Counter(data)\n',Counter(img_array)) # 调用Counter函数
bilv = []
for i in np.unique(img_array):
#print(bilv[-1]) # 对照unique数组，依次统计每个元素出现的次数
bilv.append(np.sum(img_array == i))
print(bilv[-1])
print("比率",bilv[-1]/bilv[0])
print('==========')
# 分波段读取了
if outbandsize>2:
band1 = in_ds.GetRasterBand(1) # 注意这里从1开始不是从0开始
band1_data = band1.ReadAsArray() #获取数据
band2 = in_ds.GetRasterBand(2) # 注意这里从1开始不是从0开始
band2_data = band2.ReadAsArray() #获取数据
band3 = in_ds.GetRasterBand(3) # 注意这里从1开始不是从0开始
band3_data = band3.ReadAsArray() #获取数据
compat = np.array([band1_data,band2_data,band3_data]) # compat.shape = (3, 7134, 7029)
np_remote_img = np.swapaxes(compat,0,1)
np_remote_img = np.swapaxes(np_remote_img,1,2) # np_remote_img.shape = (7134, 7029,3)
# 3 最大值最小值和唯一值
print("波段1最大值：", np.max(band1_data))
print("波段1最小值", np.min(band1_data))
print("波段1唯一值", np.unique(band1_data))
print("波段2最大值：", np.max(band2_data))
print("波段2最小值", np.min(band2_data))
print("波段2唯一值", np.unique(band2_data))
print("波段3最大值：", np.max(band3_data))
print("波段3最小值", np.min(band3_data))
print("波段3唯一值", np.unique(band3_data))
print("前3波段最大值：", np.max(np_remote_img))
print("前3波段最小值", np.min(np_remote_img))
print("前3波段唯一值", np.unique(np_remote_img))
# 读取普通图像
def read_img():
input_images_path = r"G:\微信图片.jpg"
im = Image.open(input_images_path)
array = np.array(im)
print(array.shape)
# 显示图像
plt.imshow(np_remote_img) # 显示图片
plt.axis('off') # 不显示坐标轴
plt.show()
1.2 读取文件（3种方式）
1.2.1 读取整景影像
# 将16位遥感影像转换为8位，并进行百分比截断
# https://blog.youkuaiyun.com/qq_43814272/article/details/109741119?utm_source=app&app_version=4.17.2&code=app_1562916241&uLinkId=usr1mkqgl919blen
import gdal
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from collections import Counter
import time
def read_img(input_file):
"""
读取影像
Inputs:
input_file:16位图像数据的路径
Outputs:
array_data：保存图像数据的numpy.array数组
rows:高度
cols:宽度
bands:深度
"""
in_ds = gdal.Open(input_file)
if in_ds == None:
print(input_file + "文件无法打开")
rows = in_ds.RasterYSize # 获取数据高度
cols = in_ds.RasterXSize # 获取数据宽度
bands = in_ds.RasterCount # 获取数据波段数
# 这个数据类型没搞清楚，1，2，3代表什么https://blog.youkuaiyun.com/u010579736/article/details/84594742
# GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6
datatype = in_ds.GetRasterBand(1).DataType
print("数据类型：", datatype)
# 读取方式1，全部读取
array_data = in_ds.ReadAsArray() # 将数据写成数组，读取全部数据，numpy数组,array_data.shape (4, 36786, 37239) ,波段，行，列
# img_array = in_ds.ReadAsArray(0, 0, width, height) # 将数据写成数组，分块读取，指定区域
# 读取方式三，分块读取
band1 = in_ds.GetRasterBand(1) # 注意这里从1开始不是从0开始
band1_data = band1.ReadAsArray() #获取数据
band2 = in_ds.GetRasterBand(2) # 注意这里从1开始不是从0开始
band2_data = band2.ReadAsArray() #获取数据
band3 = in_ds.GetRasterBand(3) # 注意这里从1开始不是从0开始
band3_data = band3.ReadAsArray() #获取数据
del in_ds
return array_data, rows, cols, bands
1.3 写文件（数组写为tif）
def write_img(read_path, img_array):
"""
read_path:原始文件路径
img_array：numpy数组(4, 36786, 37239) ,波段，行（高），列（宽）
"""
read_pre_dataset = gdal.Open(read_path)
if read_pre_dataset == None:
print(read_pre_dataset + "文件无法打开")
img_transf = read_pre_dataset.GetGeoTransform() # 仿射矩阵
img_proj = read_pre_dataset.GetProjection() # 地图投影信息
print("1读取数组形状", img_array.shape,img_array.dtype.name)
# GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,
if 'uint8' in img_array.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in img_array.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(img_array.shape) == 3:
img_bands, im_height, im_width = img_array.shape
else:
img_bands, (im_height, im_width) = 1, img_array.shape
filename = read_path[:-4] + '_unit8' + ".tif"
driver = gdal.GetDriverByName("GTiff") # 创建文件驱动
# 注意这里先写宽再写高，对应数组的高和宽，这里不对应才对
# https://vimsky.com/examples/detail/python-method-gdal.GetDriverByName.html
dataset = driver.Create(filename, im_width, im_height, img_bands, datatype)
dataset.SetGeoTransform(img_transf) # 写入仿射变换参数
dataset.SetProjection(img_proj) # 写入投影
# 写入影像数据
if img_bands == 1:
dataset.GetRasterBand(1).WriteArray(img_array)
else:
for i in range(img_bands):
dataset.GetRasterBand(i + 1).WriteArray(img_array[i])
2 功能大全
import gdal
import numpy as np
# 一、 读取数据
# 1 读取无人机
# 4波段
#r"J:\4影像数据大全\1无人机数据\巴南区无人机和高分影像\1巴南区无人机影像\正射影像\77-04.tif 就能读取
#r"J:\4影像数据大全\1无人机数据\巴南区无人机和高分影像\1巴南区无人机影像\正射影像\75-01.tif 不能读取
# 3波段
# H:\01HTutorWork\3GF2\3Code\case\anli12_fx\fx163_1.tif
in_ds = gdal.Open(r"H:\01HTutorWork\3GF2\3Code\case\anli12_fx\fx163_1.tif")
# 2 读取高分2号
#in_ds = gdal.Open(r"J:\4影像数据大全\3高分2号数据\2辽宁高分数据\2辽宁处理好的20景影像20210527\GF2_PMS1_E124.0_N41.9_20180924_L1A0003483014-MSS1.tif")
# 3 读取北京2号
#in_ds = gdal.Open(r"J:\4影像数据大全\2北京2号数据\1北京2号数据46景\8个县数据成果46景\TRIPLESAT_1_PSH_L4_20200415023916_002A77VI_016.pix")
# 4 读取普通的3波段jpg
in_ds = gdal.Open(r"H:\01HTutorWork\3GF2\3Code\case\anli1_feixi\2图像转为jpg\c1.jpg")
# 二、 参看数据信息，获取numpy值
print("open tif file succeed")
width = in_ds.RasterXSize # 获取数据宽度
print("宽_列：",width)
height = in_ds.RasterYSize # 获取数据高度
print("高_行：",height)
outbandsize = in_ds.RasterCount # 获取数据波段数
print("波段数：",outbandsize)
im_geotrans = in_ds.GetGeoTransform() # 获取仿射矩阵信息
print("仿射信息：",im_geotrans)
im_proj = in_ds.GetProjection() # 获取投影信息
print("投影信息：",im_proj)
datatype = in_ds.GetRasterBand(1).DataType # 这个数据类型没搞清楚，1，2，3代表什么
print("数据类型：",datatype)
# 2 分块或者分波段得到numpy数组
# img_array = in_ds.ReadAsArray(0, 0, width, height) # 将数据写成数组，分块读取，指定区域
img_array = in_ds.ReadAsArray() # 将数据写成数组，读取全部数据，结果是numpy数组
print("数组类型名",img_array.dtype.name)
print("数组形状",img_array.shape)
# 三、波段操作
# 获取特定波段的数据，参数是波段，bandnum是从1开始计数的。
dt=in_ds.GetRasterBand(3)
# 这个波段的数据
data = dt.ReadAsArray() #获取数据
data.shape
# 获取前3个波段，两种思路，一种是在gdal中拼接，生成新的栅格，一种是在numpy中转换为3维numpy数组
"""
gdal中操作参考https://blog.youkuaiyun.com/vonuo/article/details/
74783291?utm_medium=distribute.pc_relevant.none-task-blog-2~default~BlogCommendFromBaidu~default-3.essearch_pc_relevant&
depth_1-utm_source=distribute.pc_relevant.none-task-blog-2~default~BlogCommendFromBaidu~default-3.essearch_pc_relevant
numpy中操作
"""
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset # 清除内存
# 四、 数据类型转换
# 两种思路，一种是在numpy中转换，一种是在gdal中转换
"""
参考https://blog.youkuaiyun.com/vonuo/article/details/
74783291?utm_medium=distribute.pc_relevant.none-task-blog-2~default~BlogCommendFromBaidu~default-3.essearch_pc_relevant&
depth_1-utm_source=distribute.pc_relevant.none-task-blog-2~default~BlogCommendFromBaidu~default-3.essearch_pc_relevant
"""
# gdal数据类型包括
# gdal.GDT_Byte,
# gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
# gdal.GDT_Float32, gdal.GDT_Float64
#判断栅格数据的数据类型 im_data是numpy数组
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_Float32
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
# 五、 普通1或3波段图像数据的读取 它不能读取4波段数据
from PIL import Image
import numpy as np
input_images_path = r"H:\01HTutorWork\3GF2\3Code\case\anli12_fx\fx163_1.tif"
im = Image.open(input_images_path)
array = np.array(im)
array.shape
# 六、 保存
from osgeo import gdal
import os
class GRID:
#读图像文件
def read_img(self,filename):
dataset=gdal.Open(filename) #打开文件
im_width = dataset.RasterXSize #栅格矩阵的列数
im_height = dataset.RasterYSize #栅格矩阵的行数
im_geotrans = dataset.GetGeoTransform() #仿射矩阵
im_proj = dataset.GetProjection() #地图投影信息
im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵
del dataset
return im_proj,im_geotrans,im_data
#写文件，以写成tif为例
def write_img(self,filename,im_proj,im_geotrans,im_data):
#gdal数据类型包括
#gdal.GDT_Byte,
#gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
#gdal.GDT_Float32, gdal.GDT_Float64
#判断栅格数据的数据类型
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_Float32
#判读数组维数
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
#创建文件
driver = gdal.GetDriverByName("GTiff") #数据类型必须有，因为要计算需要多大内存空间
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
dataset.SetProjection(im_proj) #写入投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
if __name__ == "__main__":
os.chdir(r'D:\Python_Practice') #切换路径到待处理图像所在文件夹
run = GRID()
proj,geotrans,data = run.read_img('LC81230402013164LGN00.tif') #读数据
print proj
print geotrans
print data
print data.shape
run.write_img('LC81230402013164LGN00_Rewrite.tif',proj,geotrans,data) #写数据
import numpy as np
data = data.astype(np.float)
ndvi = (data[3]-data[2])/(data[3]+data[2]) #3为近红外波段；2为红波段
run.write_img('LC81230402013164LGN00_ndvi.tif',proj,geotrans,ndvi) #写为ndvi图像
# 七、清除内存
del dataset # 清除内存
2 参考
gdal实现对遥感影像的读写和波段选择
https://blog.youkuaiyun.com/nominior/article/details/105512346
python gdal ReadAsArray 按块读取栅格
https://blog.youkuaiyun.com/niuxuerui11/article/details/112003613
计算ndvi
https://blog.youkuaiyun.com/weixin_36080474/article/details/113538235
3 gdal包安装
https://blog.youkuaiyun.com/nima1994/article/details/79207805
扫一扫
1112
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
